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’’This  report  presents  the  results  of  an  investigation  into  the  formulation 
and  application  of  assumed- stress  hybrid  finite-elements  for  bending  of  multi- 
layer laminated  plates  and  shells.  Two  families  of  hybrid-stress-based  multi- 
layer plate  elements  are  considered;  these  elements  are  denoted  as  thick  plate 
and  moderately- thick  plate  elements.  In  the  development  of  the  thick  plate 
elements,  transverse  shear  deformation  effects  are  included  by  allowing  lines 
normal  to  the  plate  midsurface  in  the  undeformed  state  to  be  piecewise  linear 
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20  . Abstract  (Continued) 


are  included  in  an  average  sense  in  the  moderately- thick  plate  element  family  by 
assuming  that  straight  lines  normal  to  the  plate  midsurface  prior  to  deformation 
remain  straight  but  not  necessarily  normal  to  the  plate  midsurface  after  defor- 
mation. Comparison  of  the  results  obtained  by  using  the  thick  plate  and 
moderately- thick  plate  elements  with  independent  analytical  results  shows  that 
the  moderately- thick  plate  elements  are  more  efficient  and  practical  for  the 
analysis  of  multilayer  structures  having  a large  number  of  elastically  dissimila: 
layers.  » 


?oV  tiynamic  analyses,  the  Modal  Superposition  Method  (MSM)  is  employed 
to  obtain  the  timewise  solution,  and  the  Subspace  Iteration  Method  (SIM)  is 
adopted  as  an  efficient  scheme  for  calculation  of  the  lowest  few  eigenvalues 
and  eigenvectors  of  the  assembled  structure.  Both  the  SIM  and  MSM  are  programme' 
as  modules  to  be  compatible  with  a general  modular  finite-element  computer  code 
for  static  and  dynamic  analysis.  The  advantages  of  this  modular  approach  are 
demonstrated  in  a series  of  static  and  dynamic  applications  analyses. 
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SECTION  1 
INTRODUCTION 

This  report  summarizes  the  results  of  an  investigation  of  the  assumed- 
stress  hybrid  finite-element  method  [1]*  for  applications  involving  bending- 
stretching behavior  of  multilayer  laminated  composite  plates  and  shells. 

The  present  investigation  took  as  its  starting  point  the  results  of  a pre- 
vious investigation  in  which  a quadrilateral  multilayer  thick-plate  element 
was  developed  and  incorporated  in  a specialized  computer  code  for  analysis 
of  flat  plates  and  shell  segments  [2,3].  These  results  were  extended  to 
the  case  of  a triangular  element,  and  the  possibility  of  including  special- 
purpose  assumed-stress  distributions  to  satisfy  traction-free  conditions 
along  one  element  edge  was  examined  [4]. 

The  thick-plate  elements  were  found  to  be  computationally  inefficient, 
and,  hence,  were  replaced  with  moderately- thick  elements  which  were  developed 
by  extension  of  a concept  from  an  earlier  investigation  [5] . The  data 
reviewed  in  this  report  demonstrate  that  the  new  elements  are  both  computa- 
tionally efficient  and  as  accurate  as  the  thick-plate  elements  for  analyses 
of  thin  and  moderately-thick  structures,  i.e.  for  plates  with  span/thickness 
ratios  greater  than  10  and  for  shells  with  similar  ratios  of  radius-of- 
curvature  to  thickness.  These  elements  were  programmed  to  be  compatible 
with  a general  modular  finite-element  analysis  code  [6] . 

A review  of  the  material  property  matrices  and  their  transformations 
required  by  the  new  element  is  given  in  Section  2.  Sections  3 and  4 review 
the  formulations  of  the  thick  and  moderately-thick  elements,  respectively. 

The  performances  of  both  types  of  elements  in  static  stress  analyses  are 
compared  with  independent  analytical  solutions  in  Section  5. 

The  investigation  of  the  new  elements  was  also  extended  to  a formula- 
tion for  an  integrally-stiffened  quadrilateral  element,  in  which  the 
stiffener  is  treated  simply  as  an  extra  layer  which  covers  only  a part  of 
the  element's  surface.  Presented  in  Section  6 is  a comparative  evaluation 
of  this  concept  with  the  traditional  approach  in  which  stiffeners  are 
modelled  separately  as  assumed-displacement  beam-theory  elements. 


* 


Numbers  in  square  brackets  [ ] denote  references  listed  at  the  end  of  the 


Both  lumped  and  consistent  mass  matrices  were  formulated  for  the  new 


elements  (Section  4)  and  were  tested  by  means  of  dynamic  analysis  procedures 
which  were  adopted  from  the  previous  codes  [3] . The  subspace  iteration 
method  for  eigenvalue  analysis  and  the  modal  superposition  method  for 
transient  response  analysis  were  selected  on  the  basis  of  computational 
efficiency.  These  procedures  were  reprogrammed  to  be  compatible  with  a 
general  modular  finite-element  analysis  code  [6] . The  original  formulation 
of  the  modal  superposition  method  permitted  only  analyses  of  undamped 
transient  response  [3],  but  has  been  extended  in  the  present  investigation 
to  allow  damped  response  with  prescribed  modal  damping  factors  less  than 
critical.  Verification  tests  of  the  dynamic  analysis  procedures,  together 
with  some  dynamic  performance  tests  of  the  new  elements,  are  presented  in 
Section  7. 

One  phase  of  the  present  research  effort  was  the  development  of  a 
modular  dynamic  analysis  finite-element-based  computer  code.  This  has  been 
accomplished  by  modification  and  extension  of  an  existing  modular  static 
analysis  finite-element  program  [6] . A modular  code  is  arranged  as  a library 
of  subroutines,  each  of  which  performs  a series  of  computations  which,  taken 
together,  form  a single  logical  step  in  the  context  of  the  analysis.  The 
user  responsibility  now  shifts  from  the  routine  input  preparation  associated 
with  "black-box”  codes  to  the  writing  of  a MAIN  controlling  program  which 
generates  the  desired  geometrical  configuration  and  sequences  the  required 
analysis  steps.  The  use  of  a modular  approach  thus  requires  additional 
user-interaction,  but  in  turn,  provides  the  user  with  increased  flexibility 
in  comparison  with  the  use  of  "black-box"  codes. 

Incorporation  of  the  dynamic  analysis  capability  in  the  existing 
modular  code  required  that  some  modifications  be  made  in  the  original 
subroutines.  A fully  updated  user's  guide  will  be  published  in  the  near 
future  to  document  the  modified  code,  which  has  been  designated  as  FEABL-5. 
Some  preliminary  documentation  is  presented  in  Section  8 of  this  report, 
which  describes  in  detail  four  applications  programs.  The  static  applica- 
tions programs  combine  the  new  elements  with  the  existing  FEABL-2  software 
for  automatically  generated  analyses  of  a plate  with  a circular  hole  sub- 
jected to  four- point  bending,  and  stiffened  and  unstiffened  laminated 
conical  shells.  One  dynamic  application,  vibration  of  a simply- supported 
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SECTION  2 

DESCRIPTION  OF  MATERIAL  ELASTIC  CONSTANTS 


2,1  Introduction 

Each  ply  of  a composite  laminate  may  be  treated  as  a homogeneous 
orthotropic  material  for  the  purpose  of  stress  analysis.  With  homogeneity 
assumed,  there  are  at  most  9 independent  elastic  constants  for  a general 
orthotropic  material.  This  number  may  be  reduced  to  6,  of  which  5 constants 
are  absolutely  necessary,  for  the  description  of  bending-stretching  behavior 
in  thin  and  moderately  thick  laminates  inclu’ing  transverse  shear  deformation. 
Also,  two  methods  of  description  are  possible.  One,  based  on  a tensor  formula- 
tion of  the  equations  of  elasticity,  is  convenient  for  the  derivation  of  axis- 
rotation  transformations  necessary  to  translate  the  material  properties  to 
general  elastic  matrices  associated  with  arbitrarily  oriented  reference  axes. 
The  second  is  a matrix  formulation  which  is  convenient  for  the  computer- 
programming of  finite-element  stiffness  matrix  computations,  and  which  also 
corresponds  to  the  data  usually  found  in  materials  handbooks.  The  purpose 
of  this  section  is  to  review  the  relationship  between  the  two  methods  of 
description  and  their  relations  to  the  finite-element  formulations  discussed 
in  Sections  3 and  4. 

;.2  Tensor  and  Engineering  Stress  and  Strain 

The  key  to  the  relation  between  elastic  constants  is  the  relationship 
which  exists  between  the  tensor  and  "engineering"  components  of  stress  and 
strain.  The  two  systems  are  physically  identical  for  stress  (within  the 
theory  of  linear  elasticity) , the  only  difference  being  the  various  notations 
which  have  been  adopted.  The  stress  tensor  is  represented  by: 


jri 

u 

f 

r 

2 

Equivalent  "engineering"  or  vector  descriptions  of  stress  are  given  by: 


cr  = 

°k  r3  °3l  z } " 

k 

cr3  c r4 

* t iz  ^73  T21 

J 

The  braces  { } denote  a column  vector,  here  arranged  horizontally  to  save 
space . 
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The  notations  in  Eq.  2 are  commonly  used  in  the  literature  on  composite 
materials. 

In  a similar  manner,  the  strain  tensor  is  represented  by: 


% 

Kz 

A 

Xz 

Xz 

/ 

■V 

(3) 

where  the  components  y^_.  satisfy  the  continuum  strain-displacement  relations: 


X,  = ' (dUi/dxj  * Juj/foi) 


and  where  u^=u^  (x^.x^x^)  is  the  component  of  the  continuum  displacement 

field  parallel  to  axis  x..  However,  the  "engineering"  strains  z are  defined 

1 ij 

in  a slightly  different  manner: 

£,  = ^ '"l  ^31  *12  } _ ^ ■"  1 ' ? ^11  ^73  ^^23  ^ ^ (5) 

with  corresponding  strain-displacement  relations.- 

£j ■-  9ut-/JXj  +dUj-/dl{  (6) 

The  difference  arises  from  a convenience  in  the  traditional  description  of 
shear  stress-strain  relations  for  isotropic  materials,  i.e.  a.  ,*=G£. . instead 
of  T ^ =2Gy^  , where  G=E/[2(1+V)]  is  the  material  shear  modulus.  This  conven- 
tion has  been  carried  over  to  the  description  of  the  orthotropic  properties 
of  a composite  laminate,  as  will  be  seen  in  Subsection  2.4. 

2.3  Transformation  of  Stress  and  Strain 

Consider  now  two  Cartesian  axis  systems  x^x^x^  and  xix  >*3  with  different 
spatial  orientations.  Let  the  direction  cosines  between  the  two  systems  be 
given  by: 


^cos( xt , Xj) 

COJI  Xj  ( x2/ 

cos(xil 

*7>v 

r«t, 

rr\  ~ 

iz. 

3 

r? 

‘-N 

f 

Co!  (Zj,  ) 

CoiCx-,)  X2) 

Ces(*2, 

Xj) 

"it 

CoS (xu  *1  ) 

cos(Xf,  Xz) 

cos(xSl 

n3i 

”32 

Then  the  transformation  between  stress  (strain)  components  in  the  x^x^x^ 
system  and  stress  (strain  components  in  the  x^x^x^  system  can  be  expressed 
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concisely  with  the  indicial  notation  and  summation  conventions  of  tensor 
analysis: 


r. 


= mTktAj£  rkt 


V 


m7k  V X 


k/ 


(8) 


Expressions  like  Eqs.  8 are  convenient  for  compact  presentation  of 
formulations,  but  a matrix  notation  is  more  practical  for  computations.  Of 
most  interest  in  the  present  study  is  the  special  case  in  which  axes  xj'^3 
coincide,  while  there  is  a rotation  between  x^x2  anc*  ^1^2  ^ 

Eqs.  8 are  applied  to  this  case  and  "engineering"  notations  are  used,  it 
is  easy  to  show  that  0^=0,,  and: 
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and  where  0 (positive  CCW)  is  the  angle  from  x^  to  x^ . Also  of  interest  is 
the  inverse  transformation  matrix: 
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Equations  11  and  12  contain  the  familiar  Mohr  circle  transformations 
for  the  components  a1'ava12  in  the  x^xn  plane,  as  well  as  the 

corresponding  information  for  transverse  shears  °23'°31  ^^2  3*  2*"  31^  " It:  *S 
particularly  important  to  note  the  factors  of  1/2  applied  to  the  "engineering” 
shear  strains  in  Eq.  10:  the  transformation  is  valid  only  for  tensor  components. 
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a fact  which  must  be  reflected  when  "engineering"  notation  is  employed. 


2.4  Stress-Strain  Relations 


The  stress-strain  relations  for  a linear  elastic  material  are  conven- 
tionally given  in  "engineering"  form,  either  in  terms  of  compliance  constants 
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or  in  terms  of  stiffness  constants: 
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A general  orthotropic  material  possess  9 independent  elastic  constants 


[71,  such  that: 


" 


However,  Eq.  16  may  be  simplified  for  a fiber-composite  ply,  which  possesses 
more  symmetries  than  does  a general  orthotropic  material.  The  additional 
symmetries  are  physical  in  nature  and  arise  from  consideration  of  the  role 
played  by  the  fibers  in  stiffening  the  surrounding  resin  material. 

The  conventions  adopted  for  the  description  of  a typical  ply  and  its 
symmetries  are  illustrated  in  Fig.  2.  The  axes  x^x^x^  are  defined  as  the 
material  axes.  Axis  x^  is  oriented  parallel  to  the  fibers,  axis  x^  occupies 
the  in-plane  transverse  orientation,  and  axis  x^  is  oriented  perpendicular 
to  the  ply.  As  shown  in  the  figure,  the  fibers  tend  to  behave  somewhat 
' he  layers  in  their  effects  upon  elastic  response  of  the  material  in  the 
xxx2  and  x^x^  planes,  but  in  the  X2X3  plane  the  fibers  tend  to  behave  as 
small  circular  inclusions.  Thus,  there  is  good  reason  to  assume  that  the 
elastic  cross-coupling  terms  related  to  the  x^x  ) and  x^x^  planes  are  identical, 
while  those  related  to  the  X2X3  plane  have  different  values,  i.e.: 

' "^21  ’ ^22  S rS  ~ (17) 


Equation  (16)  is  thus  reduced  to  6 independent  elastic  constants: 
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The  general  notation  in  Eq.  18  is  commonly  replaced  with  a notation 
which  follows  the  conventions  for  isotropic  materials  (see  Eq.  15) . Longi- 
tudinal and  transverse  Young's  moduli  are  defined  by: 

Su  = l/E  1 S2i  - 1 /El  (19) 

The  shear  moduli  (unrelated  to  E^,E  ) are  defined  by: 

Sx}4  * 1 / 5U  = 1 / C xlz  (20) 

Poisson's  ratios  are  also  defined  by  following  the  traditional  method,  for 
example: 
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Combination  of  Eqs.  18,  19,  and  21  finally  leads  to  the  following 
expression  for  the  stress-strain  relations  in  the  material  axis  system  of 
a composite  ply: 
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where 


/Ei  ” 


The  elastic  constants  in  Eq.  22  are  commonly  identified  as  indicated  below. 
An  alternate  system  of  subscript  notation  found  in  much  of  the  literature 
on  composite  materials  is  also  shown: 

= E = Longitudinal  modulus 

E„  = E„  = Transverse  modulus 
2 T 

V>,  _ = V = Major  Poisson  ratio 
12  LT 

\>2i  = Vyj  = Minor  Poisson  ratio 

v_,  = v = Transverse  Poisson  ratio 
23  T 

* G^t  = In-plane  shear  modulus 

G„..  = G = Transverse  shear  modulus 
23  T 


Note  that  and  v^.,  are  not  independent  quantities,  in  view  of  Eq.  23. 

The  ma jor  Poisson  ratio  v is  usually  reported  when  material  properties 
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are  measured.  Values  for  and  G^  are  seldom  reported  because  of  the 

difficulty  in  measuring  these  quantities,  and  because  stress  analyses  of 
composite  laminates  have  focussed  in  the  past  on  stretching  behavior  or  on 
bending  of  thin  plates  for  which  transverse  shear  effects  can  safely  be 
neglected.  However,  reasonable  estimates  for  v and  G^  can  be  made  by 
assuming  that  these  properties  in  the  laminate  are  close  to  unreinforced 
resin  properties,  i.e.: 

V23  = + (24) 

where 


V = Resin  Poisson  ratio 
R 

E = Resin  Young's  modulus 
R 


and  where  the  resin  is  assumed  to  be  isotropic.  Some  handbook  data  for 

unreinforced  resin  properties  are  available.  Generally,  v =0.35  and 

E =10^  psi  for  many  epoxy  resins  currently  used  in  composite  materials. 
R 

If  the  compliance  matrix  S in  Eq.  22  is  inverted,  the  resulting 
stiffness  matrix  may  be  expressed  as: 
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Eq.  25  may  be  compared  with  the  corresponding  expression  for  isotropic 
material  to  further  illustrate  the  similarity  of  the  notation: 
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2.5  Stress-Strain  Relations  for  Moderately-Thick  Plates  and  Shells 


Moderately  thick  plates  and  shells  subjected  to  combined  bending- 
stretching  loads  may  be  analyzed  in  terms  of  5 components  of  stress.  The 
sixth  component  (the  stress  03  normal  to  the  plate  or  shell  midsurface)  is 
at  most  of  the  order  of  the  applied  pressure  loading,  and  can  be  neglected 
in  comparison  with  the  in-plane  components  0^,0  an<^  the  transverse  shear 

stresses  T^e  comPt-*-ance  anc*  stiffness  elastic  matrices  for  a fiber- 

composite  ply  (Eqs.  22  and  25)  then  reduce  to  the  following  forms: 
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The  compliance  matrix  S in  Eq.  27  is  in 
computation  of  assumed-stress  hybrid  element 
3 and  4) . Also,  the  in-plane  portion  of  the 
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the  form  which  is  used  for  the 
stiffness  matrices  (Sections 
elastic  stiffnesses: 
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(30) 


is  used  to  compute  assumed  stress  field  information  in  the  elements  which 

are  discussed  in  Section  4.  Note  that  neither  S nor  C require  the 

-IP 

transverse  Poisson  ratio  v„„. 
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2.6  Transformation  of  Elastic  Constants 

The  computation  of  a finite-element  stiffness  matrix  requires  application 
of  rotation  transformations  to  the  ply  elastic  constants  S and  C,  as  indicated 
in  Fig.  3.  The  figure  illustrates  a typical  ply  whose  material  reference  axes 
are  in  general  rotated  with  respect  to  the  global  reference  axes  xy.  The 
properties  S and  C are  associated  with  stress  and  strain  components  a and  e 
referred  to  the  material  axes  x^x^.  However,  it  is  necessary  to  describe  the 
behavior  of  the  finite  element  in  terms  of  the  global  stress  and  strain 
components : 

^ ^y  ^ k f £y  ^yi  ^tx  (31) 

Therefore,  global  elastic  constants: 
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must  be  obtained.  The  additional  nonzero  terms  (subscripts  61,62,54)  will 

appear  as  a result  of  the  transformation. 

G G 

Expressions  for  S and  C may  be  derived  by  first  recasting  the  matrices 
S and  C to  give  the  stress-strain  relations  in  terms  of  tensor  components 
and  then  applying  the  Mohr  circle  transformations  presented  in  Subsection  2.3. 
For  example,  e=SO  is  replaced  by: 
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where  A and  S are  related  by: 


5 - 


y 

y 

Si, 

S:z 

0 

0 

0 

Si, 

5U 

0 

0 

0 

Szi 

5zz 

0 

0 

0 

S« 

->zz 

0 

0 

0 

0 

0 

S44 

0 

0 

s*  = 

0 

0 

0 

(9 

0 

0 

% 

0 

0 

0 

0 

S55 

U o 

^0 

0 

0 

0 

0 

(9 

0 

0 

5«/2 

(34) 


In  a similar  manner. 
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while  in  the  global  reference  frame: 
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Note  that  the  matrices  S and  C are  not  symmetric.  However,  these  matrices 

~ G 

are  merely  used  to  compute  appropriate  values  for  the  symmetric  matrices  S 

Q ~ 

and  C from  which  the  final  calculations  are  made  to  obtain  the  element 


stiffness  matrix. 

Comparison  of  Fig.  3 with  Fig.  1 now  shows  that  the  Mohr  circle  trans- 
formations which  were  presented  in  subsection  2.3  apply  to  the  present  case 
G G 

with  a , y playing  the  roles  of  a,Y  in  Eqs.  9 and  10;  i.e.. 


are  given  by  Eqs.  11  and  12,  respectively.  Substitution  of 
i Eqs.  33  and  35  now  leads  to: 


Substitution  of  Eqs.  27  and  28  for  the  components  of  S*  and  C*  and  of  Eqs. 

11  and  12  into  Eqs.  40  then  leads  to  the  following  expressions  for  the 

G G 

components  of  S and  C , after  the  matrix  multiplications  have  been  carried 
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L ■ Only  the  in-plane  elastic  constants  C have  been  given  in  Eqs.  42,  since 

the  transverse  shear  stiffnesses  are  not  required  in  the  finite-element 
formulation. 
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One  phase  of  the  present  research  effort  involves  the  adaptation  of  the 
qeneral  quadrilateral  multilayer  plate  element  developed  by  Mau  [21.  In 
his  formulation,  Mau  obtained  a model  which  includes  transverse  shear  effects 
and  which  is  capable  of  representing  the  severe  cross-sectional  warping 
effects  often  observed  in  relatively  thick  laminated  plates.  This  is 
accomplished  by  assuming  that  straight  lines  normal  to  the  plate  midsurface 
prior  to  deformation  need  not  be  straight  or  normal  to  the  midsurface  of 
the  plate  following  deformation,  but  that  the  post-deformation  length  of  a 
normal  line  is  the  same  as  its  pre-deformation  length.  Mau  has  shown  that 
the  resulting  element  model  gives  accurate  displacement  and  stress  predictions 
even  for  relatively  thick  laminated  plates  where  severe  cross-sectional  warp- 


3,1  Introduction 

An  essential  ingredient  in  the  development  of  an  effective  multilayer 
plate  transient  analysis  capability  is  the  choice  of  an  accurate  and  efficient 
finite-element  model.  In  general,  one  may  choose  from  the  assumed-displacement 
model,  assumed-stress  model,  assumed-stress  hybrid  model,  or  other  mixed 
models,  but  the  choice  must  be  governed  by  the  ease  of  application,  and  the 
ability  of  the  model  to  represent  accurately  the  type(s)  of  behavior  charac- 
teristic of  the  structure  and  material  being  considered.  Earlier  studies 
by  Spilker  [5]  and  Mau  and  Witmer  [2]  suggest  that  the  assumed-stress  hybrid  model 
is  the  appropriate  choice  for  multilayer  composite  plate  structures,  and  this 
model  has  been  chosen  for  the  present  effort.  In  particular,  two  alternate 
hybrid- stress-based  multilayer  plate  elements  will  be  considered,  one  based 
on  the  approach  suggested  by  Mau  and  Witmer  [2]  (presented  in  this  section)  and 
the  other  based  on  an  improvement  of  the  approach  suggested  by  Spilker  [5] 
(presented  in  Section  4) . 


3 . 2 Review  of  Mau's  Quadrilateral  Element 


DEVELOPMENT  OF  THICK  MULTILAYER  PLATE  ELEMENTS 


SECTION  3 


F • 


ing  is  present.  It  is  believed  that  this  approach  is  the  most  general  approach 
(with  the  exception  of  a complete  three-dimensional  analysis)  for  multilayer 
plate  analyses  and  was  thus  chosen  for  the  present  study. 
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3.2.1  Derivation  of  the  Element  Stiffness  Matrix 


The  assumed-stress  hybrid  finite-element  model  is  based  on  a 
modified  complementary  energy  principle  in  which  the  requirements  of  inter- 
element traction  compatibility  and  boundary  traction  compatibility  (mechanical 
boundary  conditions)  have  been  relaxed.  The  resulting  hybrid-stress  functional, 

Tf  , may  be  stated  in  matrix  form  as 
me 

TT„c-  ?/  5 J Z'Z  £ dV  - j xVcls  + 5 T'i.  (43) 

v„  ^ 

where 

0 = stress  vector 

S = material  properties  matrix 

T = element  boundary  traction  vector 

u = element  boundary  displacement  vector 

T = prescribed  boundary  traction  vector 

V = volume  of  the  nth  element 
n 

3 V = boundary  of  the  nth  element 
n 

S = portion  of  3V  over  which  tractions  are  prescribed. 

O n 

n 

and  the  summation  in  Eq.  43  is  over  all  elements. 

The  application  of  Eq.  43  requires  the  assumption  of  a stress 
field,  a,  in  each  element  which  satisfies  the  homogeneous  equilibrium  equations, 
and  the  assumption  of  a displacement  field,  u,  along  the  element  boundary  which 
satisfies  interelement  displacement  compatibility.  For  the  present  multilayer 
element  the  stresses  in  the  mth  layer,  a™,  are  expressed  in  terms  of  a set  of 
stress  parameters,  8™,  for  the  mth  layer,  in  the  form 

<T,V'  = P ^ <44 ) 

where  P is  a function  of  the  Cartesian  coordinates  x,y,z.  The  boundary  trac- 
tions for  the  mth  layer,  T™,  can  be  related  to  the  stress  parameters,  Bm, 
using  the  direction  cosines  of  the  boundary,  and  can  be  expressed  as 


The  displacements,  u,  along  the  boundary  of  the  element  are  then  expressed 
in  terms  of  a finite  number  of  generalized  nodal  displacement  parameters 
(degrees  of  freedom) , q,  such  that  interelement  displacement  continuity  i; 
satisfied,  in  the  form 


where  L is  a function  of  position  on  the  element  boundary.  Substituting 
Eqs.  44,  45,  and  46  into  Eq.  43  yields  the  following  expression  for  tt  : 


1 1 MUmdi 


It  should  be  noted  that  the  integrations  in  Eqs.  48b  and  48d  extend  over 
the  volume  of  the  mth  layer,  V™,  and  the  boundary  surface  of  the  mth  layer, 
3v™,  respectively,  for  the  nth  element.  The  matrix  Sm  is  the  material 
property  matrix  for  the  mth  layer  in  the  x,y,z  coordinate  system  (see 
Section  2).  In  addition,  it  should  be  noted  that  the  element  H and  G 
matrices  are  "supermatrices"  composed  of  contributions  from  each  of  the  M 
layers,  a. id  3 is  a column  vector  composed  of  the  vectors  from  each  of  the 


M layers. 


The  choice  of  independent  stress  assumptions  within  each  layer 


implies  that  stresses  on  interlayer  surfaces  are  not  compatible  in  general. 
For  present  purposes,  it  will  be  assumed  that  perfect  bonding  exists  between 
layers  and  thus  the  following  stress  compatibility  requirements  should  be 
satisfied  at  the  (m+l)st  interface  (between  layers  m and  m+1) : 


^2  . S'* 


od  » ( vv»  t 

= ^xz  ) ®VZ  , 


z--  - 1 


where  h is  the  thickness  of  the  mth  layer  and  a reference  z=0  axis  is  chosen 

at  the  midsurface  of  each  layer.  In  addition  to  the  stress  compatibility 

at  interlayer  boundaries  given  by  Eq.  49,  the  requirement  that  the  stresses 

0,0,  and  0 be  zero  on  the  top  and  bottom  surfaces  of  the  plate  will 
xz  yz  z 

also  be  imposed,  i.e. 
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Equation  50a  corresponds  to  the  bottom  of  the  plate  (i.e.  bottom  surface  of 
the  first  layer)  and  Eq.  50b  corresponds  to  the  top  of  the  plate  (i.e. 
upper  surface  of  the  Mth  layer) . Equations  49  and  50  may  be  viewed  as 
constraint  conditions  on  the  full  set  of  stress  parameters,  3,  and  these 
constraint  equations  can  be  expressed  as 

A |3  •=  O (51) 

Using  the  Lagrange  Multiplier  technique,  Eq.  51  (premultiplied  by  a vector 
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of  Lagrange  multipliers,  X)  is  substituted  into  the  hybrid  functional  it 

~ n 

(Eq.  47)  to  yield  a modified  hybrid  functional,  , in  the  form 


(52) 


At  this  point  it  is  perhaps  worthwhile  to  discuss  briefly  the 

nature  and  meaning  of  tt ' . Consider  a typical  multilayer  element  as  a region 

me 

composed  of  a number  of  subregions  (in  this  case  layers) . The  application 
of  Eq.  52  can  be  shown  [8]  to  be  equivalent  to  the  use  of  the  stress  model 
of  Fraeijs  de  Veubeke  [9]  for  each  subregion  (layer)  boundary  and  the  use 
of  the  hybrid-stress  model  for  the  boundary  of  the  whole  region  (element) . 

The  implication  is  that  stress  equilibrium  must  be  satisfied  in  each  subregion 
(layer),  but  that  stress  compatibility  along  the  boundaries  of  subregions 
is  satisfied  only  in  an  average  (or  integral)  sense. 

The  expression  for  the  element  stiffness  matrix,  k,  can  now  be 
obtained  from  Eq.  52  by  setting  the  first  variation  of  with  respect  to 
8 equal  to  zero,  yielding 

H & - G ^ + AT  b ~ ° 

from  which 

fi  = N"'6  g - M'1  X 


(53) 


(54) 


Substituting  Eq.  54  into  Eq.  51  and  solving  for  the  Lagrange  multipliers, 

X,  gives 

X = (AH'V)  ~ ? ) (55) 


The  expression  which  relates  the  stress  parameters,  8>  to  the  nodal  displace- 
ments, q,  is  then  obtained  by  substituting  Eq.  55  into  Eq.  54; 

. -i 


P - [h-1^  ■ h-'at  (ah-'a1)’  A h-'s  ] g.  . B l 


(56) 


when  Eqs.  55  and  56  are  substituted  into  Eq.  52,  the  resulting  expression 


for  71’  can  be  written  in  the  familiar  form 
me 

jrJc=  % {{  fk  I - fg  } 

where  k is  the  element  stiffness  matrix  and  is  given  by 


(57) 
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W w » V* 


= C^H-'S  - (A  1''*T) 


(58) 


and  Q is  the  element  nodal  loading  vector.  It  should  be  noted  that  the  first 
term  in  Eq.  58  is  the  same  as  that  obtained  for  the  conventional  hybrid-stress 
model,  and  the  second  term  is  present  because  of  the  introduction  of  stress 
continuity  (in  an  approximate  sense)  at  interlayer  boundaries. 

3.2.2  Element  Stress  and  Displacement  Assumptions 

In  Ref.  2,  Mau  considered  various  combinations  of  linear  and 
quadratic  interpolations  for  the  stresses  in  the  interior  of  the  element 
and  displacements  along  the  element  boundary.  Based  on  the  results  of 
several  performance  tests,  he  concluded  that  the  optimum  choice  for  both 
stress  and  displacement  is  a linear  interpolation  in  x and  y. 

Mau  used  the  following  linear  interpolation  for  stresses  (which 
satisfies  the  homogeneous  equilibrium  equations  for  the  mth  layer) : 

- p?  + aT  *r ' A~V  * * C AT  * AT*  * Af  y ^ 


v 


' P 
t o 


* fiZ  - * (/C  vD - v-*") 

o-xr-  & - z c +/*?)- 

<*iy  • f A" x + A/V  + e ( a7  + a7*  * A7  y) 


(59) 


Thus,  20  stress  parameters  are  required  for  each  layer  so  that  a laminate 

composed  of  M layers  would  require  a total  of  20M  stress  parameters.  It 

should  be  noted  that  the  zero  normal  stress,  a , is  a consequence  of  the 

z 

linear  assumption  for  0,0,  and  a ; if  quadratic  terms  in  x and  y were 

x y xy 

added,  then  a would  be  nonzero.  The  P and  R matrices  for  each  layer  are 
z ~ 

obtained  from  Eqs.  59  and  the  A matrix  for  the  laminate  can  be  obtained 

from  Eqs.  59  by  applying  the  constraint  conditions  of  Eqs.  49  and  50  (for 

O and  a since  O is  zero  everywhere) . 
xz  yz  z 

The  boundary  displacement  assumption  should  be  chosen  in  such 
a way  that  a severe  cross-sectional  warping  behavior,  such  as  that  shown  in 
Fig.  4,  can  be  approximately  represented.  This  can  be  accomplished  by  defin- 
ing the  inplane  displacement  behavior  in  terms  of  translational  nodal 
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displacement,  w,  is  assumed  to  be  linear  along  each  side  of  the  element 
but  constant  in  the  z direction  for  all  layers.  Therefore,  the  total 
number  of  degrees  of  freedom  at  a node  consists  of  one  w and  2(M+1) 
inplane  displacements  u and  v,  where  M is  the  total  number  of  layers.  The 
total  number  of  degrees-of-freedom  for  the  present  four-node  general  quadri- 
lateral element  is  thus  8M+12. 

For  a typical  layer,  m,  the  displacement  interpolation  along 
side  1-2  (between  nodes  1 and  2 — see  Fig.  5),  for  example,  is  given  by 

uv,  - { (u7"+ur)('-s)  * i («r'*ur)s 

--  i (vr*V  <")(<--'>  + ^('r~*,4.vr,~)s  [(vrr'i  V.-Xi-S1!  *UC".  (60) 

^,7  - w,  C i - s')  + s 

where  s is  a nondimensional  parameter  which  takes  on  the  values  s=0  at 
node  1 and  s=l  at  node  2.  Expressions  similar  to  Eqs.  60  can  be  obtained 
for  each  layer  and  each  side  of  the  element. 

In  practice,  the  supermatrices  H and  G are  obtained  by  forming 
the  matrices  Hm  and  G™  for  the  mth  layer  and  inserting  the  appropriate  terms 
in  A corresponding  to  the  mth  layer.  When  all  layers  have  been  processed, 
the  H,  G,  and  A matrices  are  complete.  The  element  B matrix  can  now  be 
obtained  from  Eq.  56  and  the  element  stiffness  matrix  can  be  obtained  from 
Eq.  58.  In  subsequent  discussions,  the  4-node  general  quadrilateral  multi- 
layer plate  element  based  on  Eq.  52  and  utilizing  the  stress  assumption  of 
Eq.  59  and  displacement  assumption  of  Eq.  60  will  be  termed  element  ELEMZ. 

3.2.3  Derivation  of  the  Element  Mass  Matrix 

As  was  shown  in  the  previous  subsections,  the  use  of  the  hybrid- 
stress  finite-element  model  for  static  analysis  is  equivalent  to  the  use  of 
the  conventional  assumed-displacement  model  in  the  sense  that  the  resulting 
matrix  equations  are  expressed  in  terms  of  unknown  nodal  displacement  parameters. 
The  hybrid-stress  model  thus  represents  an  alternate  way  of  obtaining  an  element 
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stiffness  matrix.  For  dynamic  analysis,  the  corresponding  element  mass 
matrix  is  required. 

The  derivation  of  the  element  mass  matrix  is  most  conveniently 
seen  by  considering  a functional  in  the  form  of  a Hellinger-Reissner 
principle  [10,11]  for  the  free  vibration  of  a continuum, 


Tf^  = f [ I £'i  ®*£j<rkJ,  +^<r;j  (cn.j 

- J T (u.-uj  <r/s  ^ 


In  Eq.  61,  tensor  notation  and  the  summation  convention  have  been  employed. 
Also,  p is  the  material  density,  is  the  displacement  field  in  the  interior 

of  the  element,  u.  is  the  displacement  field  on  the  boundary  of  the  element 

- • 

(note  that  u.  need  not  be  equal  to  u . on  3v  ),  u.  is  the  velocity  field  in 
x i n l 

the  interior  of  the  element,  and  a comma  denotes  partial  differentiation. 

The  reduction  of  7 1 to  correspond  to  the  hybrid-stress  model  is  accomplished 

mR 

by  assuming  that  the  boundary  tractions,  T.,  are  related  to  the  stress,  a.  ., 


Ti  - S Oj 


on 


where  V.  is  the  direction  cosine  tensor  on  the  element  boundary,  3v  . Then 
3 n 

the  Divergence  theorem 


1 k S - I - I S/.j  M‘ 


is  applied  to  Eq.  61  to  yield 


" § f J £"Z  S'  ^ - °*j.j  U'  - i ]dV 

r l (6  ■ 

+ ) Tiu,  as  \ 

If  the  stress  field  within  each  element  exactly  satisfies  the  homogeneous 
portion  of  the  equilibrium  equations,  i.e. 
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then  Eq.  64  can  be  rewritten  as 


? [i  J ^ ^ dv  + i dv/  - JT^  d 

V*  a>v„ 

Equation  66  may  be  viewed  as  a modified  hybrid-stress  functional  so  that  the 
governing  functional  for  dynamic  analysis  can  be  written  in  matrix  form  as 


7/~ 


(67) 


, - ? ^-7-/crTS(r<rly*7^p^T^civ/_r-rT5c/s  -»  T t Ta  cli  l 
" v-  ^ ~ J 

Equation  67  is  the  dynamic  equivalent  of  Eq.  43,  and  it  should  be  clear 
that  the  second  term  in  Eq.  67  will  yield  the  desired  element  mass  matrix. 
As  is  done  for  static  applications,  the  stresses  are  expressed  in  terms  of 
unknown  stress  parameters,  (3,  and  the  boundary  displacements,  u,  are 
expressed  in  terms  of  unknown  nodal  displacement  parameters,  q.  In 
addition,  a displacement  assumption  for  u in  the  interior  of  the  element 
is  now  required  in  terms  of  the  nodal  displacement  parameters,  q, 


U - N % 


(68) 


It  should  be  recalled  that  the  interior  displacement  assumption  need  not  be 

equal  to  the  boundary  displacement  assumption  (i.e.  u need  not  be  equal  to 

u on  9v^) . Then  Eqs.  44,  45,  46,  and  68  are  substituted  into  Eq.  67,  and 

the  constraint  condition  (Eq.  51)  on  stresses  at  interlayer  boundaries  is 

introduced.  Finally,  the  stress  parameters,  3«  and  Lagrange  multiplier,  X, 

are  eliminated  in  favor  of  the  nodal  displacement  parameters,  q,  the  resulting 

expression  for  tt ’ being  qiven  by 
me 


- 2l  l { fk  I * i i - j 


(69) 


where  k is  the  element  stiffness  matrix  given  by  Eq.  58  and  m is  the  element 
mass  matrix  given  by 

vvi  - \ P ^ ~ ^ (70) 

V* 

It  should  be  noted  that  the  use  of  Eq.  67  is  not  fully  consistent  with  the 
hybrid-stress  model.  This  is  because  the  hybrid-stress  model  requires  that 
the  equilibrium  equations  be  satisfied.  However,  for  dynamic  analysis,  the 
equilibrium  equations  include  an  inhomogeneous  portion  corresponding  to 


inertia  terms,  and  are  thus  not  satisfied  by  the  stress  assumptions  in  the 
present  formulation.  Since  only  the  homogeneous  equilibrium  equations  are 
satisfied,  it  is  more  appropriate  to  view  Eq.  67  as  a modified  Hellinger- 
Reissner  functional  or  a modified  hybrid-stress  functional.  For  convenience, 
the  element  mass  matrix  derived  from  Eq.  67  will  be  termed  a "hybrid  rational" 
mass  matrix. 

Now  consider  the  development  of  a hybrid-rational  mass  matrix 
corresponding  to  element  ELEMZ.  Since  the  element  boundary  displacements 
for  ELEMZ  are  assumed  to  be  linear,  a convenient  and  suitable  interpolation 
for  N would  be  a bilinear  expansion  in  terms  of  a pair  of  transformed  coordinates 

(S.n) : 


\T 


w 


<rf\ 


i { [i  (<”'<•  -c)  ♦ ''■rK } 

l-l  J 

^2  N,  | 

i=i  -* 


(71) 


where 

* (/-  JY  i - 7 ) 
a/2  - ( l-  'i ) f 


Equations  71  have  been  written  for  the  mth  layer  of  the  multilayer  element, 
and  it  should  be  noted  that  the  original  coordinates  (x,y)  are  related  to 
the  transformed  coordinates  (£,n)  by 
4- 


* - 1 

tM  (73) 

4 

v - 2.  Vi 

t=  I 


where  (x.,y.)  are  the  coordinates  of  the  ith  node  of  the  element.  The  element 
i i 

mass  matrix  is  obtained  from  Eq.  70  by  integration  over  the  volume  of  each 


layer  of  the  laminate,  with  N and  p in  Eq.  70  replaced  by  the  particular  Nm 
and  pm  corresponding  to  the  mth  layer. 

3.3  Adaptation  of  the  Mau  Element  to  Triangular  Shape 

In  many  structural  problems  of  interest,  there  are  regions  of  the  struc- 
ture in  which  high  stress  gradients  may  exist,  and  other  regions  where  stress 
gradients  are  low.  In  practice  it  is  most  advantageous  to  be  able  to  use  a 
refined  finite-element  mesh  in  regions  of  high  stress  gradients  and/or  other 
regions  of  interest,  and  a more  coarse  mesh  arrangement  in  regions  of  low 
stress  gradients.  To  do  this,  a suitable  element  must  be  available  which 
can  be  used  as  a transition  element  between  the  refined  and  coarse  mesh 
regions,  and  which  is  compatible  with  the  elements  in  these  regions.  The 
transition  element  which  will  be  employed  in  the  present  effort  is  a triangular- 
shaped multilayer  plate  element. 

Because  this  triangular  element  is  to  be  used  in  conjunction  with  the 
quadrilateral  element  ELEMZ  described  in  Subsection  3.2,  the  primary  restric- 
tion is  that  the  displacement  assumption  along  the  element  boundary  and  the 
number  of  generalized  displacement  parameters  at  a node,  for  the  triangular 
element,  must  be  identical  to  those  used  for  ELEMZ.  Thus,  the  development 
of  the  present  triangular  plate  element  is  based  on  the  same  boundary  dis- 
placement and  interior  stress  assumptions  as  used  in  ELEMZ.  Because  of  this, 
only  straightforward  modifications  of  the  volume  and  surface  integrations 
used  to  generate  ELEMZ  are  necessary  to  obtain  the  triangular  element. 

Mau  has  included  in  Ref.  2 a comparison  of  results  obtained  for  a 
single-layer  isotropic  plate  problem  by  using  several  different  element 
types  (quadrilateral  and  triangular,  with  a variety  of  different  stress  and 
displacement  assumptions)  including  single-layer  versions  of  the  present 
quadrilateral  and  triangular  multilayer  plate  elements.  Based  on  these 
initial  comparisons,  Mau  concluded  that  the  present  quadrilateral  element 
yields  reasonably  accurate  solutions,  whereas  the  present  triangular  element 
yields  results  which  are  too  stiff  for  practical  use  of  the  triangular  element 
to  model  the  entire  plate.  However,  when  used  solely  as  a transition  (mesh 
expander)  element,  the  stiffening  effects  of  the  triangular  element  should 
be  less  pronounced  in  terms  of  determining  the  overall  structural  behavior. 
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3.4  Development  of  a Quadrilateral  Element  with  a Traction-Free  Edge 

When  the  finite-element  method  is  employed  for  the  analysis  of  struc- 
tures containing  cutouts,  it  is  often  necessary  to  use  a highly  refined 
mesh  in  the  region  of  the  cutout  so  that  the  stress  gradient  near  the 
cutout  and  the  traction-free  condition  at  the  surface  of  the  cutout  may  be 
more  closely  approximated.  Pian  has  presented  some  results  for  single  layer, 
isotropic  materials  which  suggest  that  improved  results  may  be  obtained  near 
traction- free  boundaries  by  using  a special  assumed-stress  hybrid  element 
which  exactly  satisfies  the  traction-free  condition  [12].  Thus,  one  phase 
of  the  present  effort  was  the  investigation  and  evaluation  of  a quadrilateral 
shaped,  multilayer  plate  element  for  which  the  traction-free  condition  is 
exactly  satisfied  on  one  of  the  element  boundaries.  In  order  for  such  an 
element  to  be  of  practical  use  in  engineering  analyses  it  must  be  computa- 
tionally efficient  and  it  must  yield  displacement  and  stress  results  near 
a cutout  which  are  significantly  better  than  the  results  which  could  be 
obtained  by  using  conventional  elements. 

In  practice,  because  this  special  element  may  be  linked  to  ELEMZ,  the 
boundary  displacement  assumption  employed  for  the  traction-free  element 
must  be  identical  to  the  boundary  displacement  assumption  employed  in  ELEMZ. 
The  stress  assumption  chosen  for  the  traction-free  element  must  satisfy  the 
equilibrium  equations  and,  in  addition,  must  yield  zero  surface  tractions 
(not  necessarily  zero  stresses)  when  evaluated  on  the  traction-free  boundary 
of  the  element.  In  essence,  the  only  difference  between  the  formulation  of 
the  traction-free  element  and  the  formulation  of  ELEMZ  is  in  the  choice  of 
an  appropriate  stress  field  in  the  interior  of  the  element. 

Several  different  stress  assumptions  were  considered  in  a study  reported 
in  Ref.  4.  For  the  sake  of  completeness,  the  essence  of  this  study  will  be 
summarized  here.  For  convenience,  the  notation  A^A^A^A^  w-’tH  be  adopted  to 
describe  the  order  of  the  stress  assumption.  In  this  notation,  A^,  will 
denote  the  highest  order  of  z included  in  the  stress  assumption  (either 
linear,  L,  or  quadratic,  Q) . will  denote  the  highest  order  (either  L 

or  Q)  of  x and  y contained  in  the  portion  of  the  stress  assumption  of  order 

z°  (constant) . A will  denote  the  highest  order  (either  L or  Q)  of  x and  y 

3 1 
contained  in  the  portion  of  the  stress  assumption  of  order  z (linear) . If 
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A =Q,  then  A^  will  denote  the  highest  order  (either  L or  Q)  of  x and  y 

2 

contained  in  the  portion  of  the  stress  assumption  of  order  z (quadratic) . 

It  should  be  noted  that  this  notation  applies  only  to  the  inplane  stress 

assumption  (a  , a ,a  ),  and  it  is  assumed  that  the  interpolation  is  complete 
x y xy 

up  to  the  order  stated  (i.e.  L implies  that  terms  of  the  order  1,  x,  and  y 

2 2 

are  present,  and  Q implies  that  terms  of  the  order  l,x,y,xy,x  , and  y are 
present).  For  example,  the  stress  assumption  used  for  ELEMZ  (Eq.  59)  would 
be  denoted  by  LLL  (i.e.  linear  in  z,  and  linear  in  x and  y in  the  z°  and  z^ 
terms) . 

In  order  to  obtain  the  required  stress  assumption,  an  assumption  is 

made  for  the  inplane  stress  components  (0  ,a  ,a  ) and  the  remaining  stress 

x y xy 

components  (O^  2'az)  are  obtained  by  solving  the  equilibrium  equations. 

The  resulting  stress  assumption  thus  satisfies  the  homogeneous  equilibrium 
equations  as  required.  Next  the  traction-free  conditions  must  be  satisfied. 

For  convenience,  a local  (element)  axis  system  is  adopted  in  such  a way  that 
the  traction-free  edge  coincides  with  the  x axis;  this  simplifies  the  traction- 
free  conditions  to 
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(74) 
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The  conditions  given  by  Eqs.  74  are  then  imposed  on  the  stress  distribution 

and  appropriate  stress  parameters  (B's)  are  eliminated  so  that  Eqs.  74  are 

exactly  satisfied.  The  resulting  stress  assumption  has  a reduced  number  of 

stress  parameters  and  satisfies  both  equilibrium  and  the  traction-free 

condition  exactly.  It  should  be  noted  that  the  resulting  number  of  B's 

following  reduction  cannot  be  arbitrarily  small.  Tong  and  Pian  [13]  have 

shown  that  the  number  of  stress  parameters,  n0,  must  be  greater  than  or  equal 

p 

to  the  number  of  nodal  displacement  parameters,  n^,  minus  the  number  of  rigid 

body  modes,  n (=6  for  the  present  elements),  i.e. 

R 

>'  (?5) 
Equation  75  is  a necessary  but  not  sufficient  condition  to  guarantee  the 
existence  of  a solution.  Even  when  Eq.  75  is  satisfied,  so-called  "additional 


kinematic  modes"  may  be  present  [5,8].  These  modes  behave  in  a fashion 
similar  to  rigid-body  modes;  if  they  are  not  properly  constrained*,  the 
assembled  stiffness  matrix  will  either  be  singular  or  so  poorly  conditioned 
that  reliable  results  cannot  be  obtained.  A more  thorough  discussion  of 
this  phenomenon  as  well  as  procedures  for  determining  these  modes,  if  present, 
may  be  found  in  Refs.  5 and  8. 

Initially,  attempts  were  made  to  choose  stress  assumptions  which  had 

a minimum  number  of  stress  parameters.  The  first  assumption  attempted,  LLL, 

coresponds  to  that  used  in  ELEMZ.  However,  after  Eqs.  74  were  satisfied, 

it  was  found  that  Eq.  75  was  not  satisfied;  thus  the  LLL  assumption  could 

not  be  used.  Next,  the  order  of  the  bending  stress  behavior  was  increased 

to  yield  a LLQ  assumption  which  contained  36  stress  parameters.  After 

application  of  Eqs.  74,  the  number  of  stress  parameters  was  reduced  to  21 

(which  satisfies  Eq.  75  even  for  a single  layer) . However,  numerical 

-1  T 

difficulties  were  encountered;  after  assembly,  the  AH  A supermatrix 
was  found  to  be  singular  and  thus  could  not  be  inverted  as  needed  to  obtain 
the  element  stiffness  matrix  (see  Eq.  58) . By  tracing  the  origin  of  this 
ill  behavior,  it  was  found  that  a linear  stretching  part  in  the  assumed 
stress  field  does  not  furnish  the  odd  powers  of  z necessary  to  make  the 
interlayer  conditions  A3=0  linearly  independent.  The  lack  of  odd-power 
terms  gives  an  identical  relation  for  the  top  and  the  bottom  of  each  layer. 

In  other  words,  the  shear  stresses  appear  to  be  the  same  at  the  bottom  and 
the  top  of  each  layer  and  from  layer  to  layer,  which  is  too  restrictive. 

Another  way  to  increase  the  number  of  $'s  is  to  use  an  in-plane  stress 
field  linear  in  x and  y as  in  ELEMZ  (computationally  a good  bargain)  but 
quadratic  in  z.  This  third  formulation  (QLLL)  satisfies  the  "odd  powers" 
condition  and  Eq.  75  for  two  or  more  layers,  but  is  found  to  contain  addi- 
tional kinematic  modes.  Here,  the  approximation  is  probably  too  crude- 

only  13  S' s,  no  a , no  O , no  a . In  addition,  increasing  the  order  in 
z y yz 

z while  keeping  the  x y behavior  linear  does  not  appear  logical,  since  the 
xy  dimensions  of  the  element  are  generally  larger  than  its  thickness.  More- 
over, Pagano's  results  [14]  for  the  cylindrical  bending  of  multilayer  plates 
show  a quasi-linear  behavior  of  the  in-plane  stresses  in  the  z direction. 


The  structure  must  be  constrained  in  such  a way  that  all  deformation  modes 
of  the  structure  corresponding  to  the  additional  kinematic  mode  are  con- 
strained. For  most  cases,  this  is  accomplished  simply  by  imposing  the 
physical  boundary  condition  of  the  problem. 


4 


29 


The  final  stress  assumption  attempted  contained  quadratic  expansions 
in  x and  y in  both  the  stretching  and  bending  portions  of  the  inplane  stress 
assumption  (i.e.  LQQ) . After  reduction,  the  stress  assumption  contained 
25  0's.  However,  the  assembled  stiffness  matrix  was  found  to  be  singular 
when  a simple  plate  analysis  was  performed.  The  source  of  this  difficulty 
(currently  under  investigation)  appears  to  be  the  presence  of  a single  addi- 
tional kinematic  mode. 

In  view  of  the  various  difficulties  encountered  in  the  development  of 
a traction-free  element,  no  further  stress  assumptions  were  attempted.  However, 
this  was  not  the  sole  reason  for  abandoning  the  search  for  a traction-free 
element.  In  particular,  it  was  found  that  the  use  of  ELEMZ  was  impractical 
for  engineering  analysis  because  of  restrictions  on  the  number  of  layers  which 
can  be  accommodated,  as  explained  in  the  next  subsection.  In  addition,  per- 
formance tests  (Section  5)  indicated  that  an  alternate  multilayer  plate  element 
(Section  4)  would  prove  to  be  more  practical  for  use  in  engineering  analysis, 
and  could  also  yield  reasonable  estimates  of  stress  distribution  near  a cutout. 

In  summary,  it  is  believed  that  a stress  assumption  of  higher  order  than 
the  LQQ  assumption  will  be  required,  which  suggests  that  the  number  of  stress 
parameters  after  reduction  will  exceed  25  B's  per  layer.  Because  of  core 
storage  requirements  associated  with  the  ELEMZ  formulation,  such  an  element 
would  be  restricted  to  structures  having  few  layers.  It  is  believed  that 
such  an  element  would  be  of  little  use  in  practical  engineering  analyses. 

3.5  Discussion 

The  elements  developed  in  this  section  are  based  on  the  assumption  of 
an  independent  set  of  stress  parameters,  0,  in  each  layer  and  nodal  displace- 
ment parameters,  q,  which,  in  effect,  allow  for  piecewise  linear  rotation  of 
the  cross-section  of  each  layer.  This  level  of  generality  may  be  necessary 
for  thick  laminated  plates  exhibiting  severe  cross-sectional  warping.  The 
purpose  of  this  subsection  is  to  discuss  the  price  to  be  paid  for  including 
these  generalities.  The  discussion  will  center  on  element  ELEMZ,  but  the 
observations  are  valid  for  the  other  elements  based  on  the  same  formulation. 

The  first  penalty  is  in  terms  of  computer  core  storage  requirements. 
Assuming  a laminate  composed  of  n layers,  the  total  number  of  stress  parameters 
required  would  be  20n  and  the  total  number  of  nodal  degrees  of  freedom  would 

be  (8n+12 ) for  a single  element  ELEMZ.  The  G supermatrix  (Eq.  48c)  would 

2 ~ 

'equire  (160n  +240n)  words  of  storage,  the  H supermatrix  (Eq.  48a)  would 
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require  (210n)  words  of  storage  (storing  only  the  lower  triangle  of  each  of 

the  symmetric  Hm  matrices),  the  element  B matrix  (Eq.  56)  would  require 
2 

(160n  +240n)  words  of  storage,  and  the  element  stiffness  matrix,  k (Eq.  58) , 

2 

would  require  (32n  +100n+72)  words  of  storage  (storing  only  the  lower  triangle 

of  the  symmetric  stiffness  matrix) . The  storage  of  all  of  these  arrays  in 

2 

core  would  require  (352n  +790n+72)  words  of  storage.  Tabulated  below  are 
the  storage  requirements  in  words  and  BYTES  (assuming  4 BYTES/word  for  single 
precision  arithmetic)  for  several  values  of  n: 

n Storage  (words)  Storage  (bytes) 


1 

1,214 

4,856 

3 

5,610 

22,440 

6 

17,484 

69,936 

9 

35,694 

142,776 

15 

91, 122 

364,488 

20 

156,672 

626,688 

30 

340,572 

1,362,288 

Thus,  based  solely  on  in-core  storage  requirements  for  the  generation  of  a 
stiffness  matrix,  certain  limitations  will  exist  on  the  maximum  number  of 
layers  which  can  be  accommodated. 

In  many  engineering  analyses,  the  largest  single  block  of  storage 
required  is  for  the  assembled  stiffness  matrix  (and  mass  matrix,  if  consistent 
mass  matrices,  rather  than  lumped  mass  matrices,  are  employed) . For  illustra- 
tive purposes,  assume  that  ELEMZ  is  to  be  used  to  analyze  a square  flat  plate 

which  is  modeled  by  a uniform  mesh  having  m elements  in  the  x-direction  and 

2 2 

m elements  in  the  y-direction  (i.e.  a total  of  m elements,  and  (m+1)  nodes). 
Assume  that  the  nodes  are  numbered  optimally  to  yield  the  smallest  bandwidth 
in  the  assembled  stiffness  matrix,  and  that  the  lower  triangle  of  the  assembled 
stiffness  matrix  is  stored  row  by  row  (in  a vector)  starting  in  each  row  with 
the  first  nonzero  term.  The  storage  requirements  can  be  estimated  as  the 
total  number  of  degrees  of  freedom  times  the  average  semi-bandwidtli . For 

ELEMZ,  each  node  has  (2n+3)  degrees-of- freedom  (n=number  of  layers) , so  that 

2 

the  total  number  of  degrees-of- freedom  would  be  (m+1)  (2n+3)  and  the  average 

semi -bandwidth  can  be  estimated  as  (m+2) (2n+3) . Thus,  the  total  number  of 
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words  of  storage  required  for  the  assembled  stiffness  matrix  is  given  approxi 

2 2 

mately  by  (m+1)  (m+2) (2n+3)  . Tabulated  below  are  the  storage  requirements 

(in  thousands  of  words)  for  several  values  of  n and  m: 


These  requirements  would  double  for  dynamic  analyses  if  both  the  assembled 
stiffness  and  mass  matrices  were  stored  in  core  simultaneously.  As  a result 
of  these  storage  requirements,  use  of  ELEMZ  may,  in  practice,  be  limited  to 
laminated  plates  composed  of  fewer  than  five  layers  if  a complete  in-core 
solution  is  desired. 

The  second  penalty  for  using  element  ELEMZ  is  in  terms  of  computation 

time.  Typically,  tl  most  time  consuming  operation  in  the  formation  of  an 

element  stiffness  matrix  by  the  hybrid-stress  method  is  the  inversion  of  the 

H matrix.  For  element  ELEMZ  the  20  by  20  Hm  matrix  (associated  with  the  mth 

layer)  must  be  inverted  for  each  of  the  n layers  in  the  laminate.  In  addition, 

-IT 

a single  inversion  of  the  (AH  A ) supermatrix  is  required.  The  computation 
time  required  to  process  the  assembled  matrix  equations  will  increase  rapidly 
as  the  total  number  of  degrees  of  freedom  increases,  which  is,  in  turn,  related 
to  the  total  number  of  layers. 

In  view  of  these  considerations,  an  alternate  element  formulation  may 
be  required  for  which  the  storage  requirements  are  less  severe,  while  main- 
taining an  acceptable  degree  of  accuracy.  Such  a formulation  is  presented 
in  the  next  section. 


SECTION  4 


DEVELOPMENT  OF  MODERATELY- THICK  MULTILAYER 
PLATE  ELEMENTS 


4.1  Introduction 

In  this  section,  an  alternate  formulation  will  be  presented  for  a multi- 
layer plate  element  based  on  the  hybrid-stress  model.  In  the  formulation  of 
element  ELEMZ  (Section  3) , the  number  of  degrees  of  freedom  per  element 
was  dependent  on  the  number  of  layers  so  that  severe  cross-sectional  warping 
effects  from  layer  to  layer  could  be  approximately  represented.  This,  in 
turn,  dictated  that  the  number  of  stress  parameters  per  element  also  be 
dependent  on  the  number  of  layers  so  that  the  B-q  relation,  discussed  in 
Subsection  3.4,  would  be  satisfied  for  an  arbitrary  number  of  layers.  As 
a result,  the  computer  storage  requirements  when  using  ELEMZ  may  become 
prohibitively  large  even  for  relatively  few  layers.  To  alleviate  these 
difficulties,  a less  general  deformation  behavior  is  incorporated  in  the 
present  alternate  element  formulation.  In  particular,  it  is  assumed  that 
straight  lines  normal  to  the  plate  midsurface  prior  to  deformation  remain 
straight  (but  not  necessarily  normal  to  the  plate  midsurface)  after  deforma- 
tion (Fig.  6).  By  utilizing  this  more  restrictive  assumption,  transverse  shear 
deformation  effects  are  still  included  (because  of  the  assumption  of  non- 
normal rotations  of  the  cross-section) , but  the  layer  to  layer  cross-sectional 
warping  will  be  modeled  only  in  an  average  sense.  It  is  believed  that  such 
an  assumption  is  adequate  for  thin  or  moderately  thick  laminates. 

Such  an  element,  based  on  the  hybrid-stress  model,  was  developed  by 
Spilker  [51.  In  the  earlier  element,  the  number  of  degrees  of  freedom  was 
fixed  at  five  per  node  (independent  of  the  number  of  layers)  and  the  number 
of  stress  parameters  for  the  laminate  was  18+2n,  where  n is  the  total  number 
of  layers.  The  calculation  of  the  element  stiffness  matrix  uses  the  same 
expression  (Eq.  58)  as  used  for  ELEMZ.  Element  performance  tests  cited  in 
Ref.  5 (displacement  results  only)  suggest  that  reliable  results  can  be 
obtained  for  thin  and  moderately- thick  plates  where  transverse  shear  deforma- 
tion effects  are  present. 

The  present  alternate  element  is  based  on  an  improvement  of  the  formula- 
tion given  in  Ref.  5.  As  will  be  shown,  the  number  of  stress  parameters  and 


33 


nodal  degrees  of  freedom  in  the  resulting  element  are  independent  of  the 
number  of  layers  in  the  laminate,  and  the  expression  for  the  element  stiff- 
ness matrix  is  identical  to  that  employed  in  the  conventional  assumed-stress 
hybrid  model. 

4.2  Derivation  of  Stiffness  and  Mass  Matrices  for  a Quadrilateral  Element 

The  derivation  of  the  stiffness  and  mass  matrices  for  the  present  four- 
node  general  quadrilateral  multilayer  plate  element  follows  the  development 
given  in  Section  3 with  the  following  key  simplifications: 

1.  The  stresses  within  each  layer  are  assumed  in  terms  of  a set  of 
stress  parameters,  3,  which  are  the  same  for  each  layer,  and  thus 
correspond  to  the  laminate  as  a whole. 

2.  In  addition  to  satisfying  the  homogeneous  equilibrium  equations, 
the  stress  assumption  must  exactly  satisfy  the  interlayer  stress 
compatibility  relations  given  by  Eqs.  49  and  50. 

3.  The  displacement  behavior  of  the  laminate  is  defined  in  terms  of 
translational  displacements,  u,v,  and  w (in  the  x,y,  and  z direc- 
tions, respectively)  of  the  plate  midsurface  and  rotations  0^  and 
9 of  lines  normal  to  the  midsurface  prior  to  deformation.  The 

y 

rotations  0^  and  9^  correspond,  but  are  not  identical  to  3w/9y  and 
9w/3x,  respectively. 

The  introduction  of  the  second  simplification  implies  that  Eq.  51  is 

exactly  satisfied  and  need  not  be  introduced  into  the  hybrid  functional. 

As  a result,  the  conventional  form  of  IT  (Eq.  43)  can  be  used  and  the 

me 

relation  between  stress  parameters,  3,  and  nodal  displacement  parameters, 
q,  (Eq.  56)  reduces  to 

/3  --  H''&  $ - g ? (76) 

The  expression  for  the  element  stiffness  matrix,  k,  (Eq.  58)  then  reduces 
to 

^ = $ ~ ~ (77) 
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Equations  76  and  77  are  the  same  as  those  obtained  in  the  conventional 
hybrid-stress  model.  For  the  present  element,  the  stresses  in  the  mth 
layer  are  expressed  in  terms  of  an  interpolation  matrix  for  the  mth  layer, 
Pm,  and  a set  of  stress  parameters,  B>  which  are  the  same  for  each  layer. 
Thus,  the  stresses,  om,  and  boundary  tractions,  Tm,  in  the  mth  layer  will 
be  assumed  in  the  form: 


f £ 

T™--  R”  £ 


The  element  H and  G matrices  are  then  defined  as 


H --  H ' + H 2 + 


^ / - Z ^ M 

Cq  -=  Qf  -f-  Gr  +•••/-<% 


5™  p-  dv 


j S-’ t 


the  integration  extending  over  the  volume,  V™,  and  boundary,  3V™,  of  the  mth 
layer  of  the  nth  element.  Because  of  the  first  and  third  simplifications, 
the  dimensions  of  these  matrices  are  fixed  (i.e.  independent  of  the  number 
of  layers) . The  method  for  determining  the  P™  matrices  will  now  be  discussed. 

It  is  required  that  the  stress  assumption  in  each  layer  be  related  to 
a set  of  stress  parameters,  3,  which  correspond  to  the  entire  laminate.  To 
accomplish  this,  the  inplane  strains  for  the  entire  laminate  are  first 
expressed  in  terms  of  a set  of  parameters,  3,  using  a linear  interpolation: 


f3,  + P+*  + fill  + Z t-  fro*  &3*  *?!<■  y ; 
pz + /V  * AjY  + 2 ( £/ + £/«*  * £7  y ) 

= + £.*  + M + * ( £2  + As  X ^ As  Y ^ 
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where  a reference  axis  z=0  defines  the  midsurface  of  the  laminate.  The 
inplane  stresses  in  the  mth  layer  are  related  to  the  inplane  strains  by 
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(83a) 


(83b) 


(83c) 


The  Aj  (i»j-l/2,6)  terms  are  the  inplane  material  stiffness  coefficients 
transformed  into  the  global  x,y,z  coordinate  system,  as  given  in  Section  2 


(Eqs.  42).  The  remaining  stress  components  for  the  mth  layer  are  obtained 
from  the  equilibrium  equations  as: 

<r  - - 5<<c  + <y)d* 

■ - * [ c"p*  A,  4 a:  a 4 Ca  As  4 C Ay  I (83d) 
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- * [ c”  p7  4 /S8  + c-  p,  + c,?  pA  + c~  ^ + cZ  pk  ] 

■ f [ cu  At  * c"  a., + c ar  + c a, 3 + c a«  + c pj + p; 


^Xj,*  + <3"-V£>'J  ) tiZ  =0 


where  D™  and  D™  are  the  constants  of  integration  for  the  mth  layer.  The 

stress  assumption  for  the  mth  layer  given  by  Eqs.  83a  through  83f  satisfies 

equilibrium  as  required.  Also,  it  should  be  noted  that  the  zero  am  is  a 

z 

consequence  of  the  linear  (in  x and  y)  assumptions  chosen  for  the  inplane 
stres  es. 

For  the  present  formulation,  the  stress  assumption  is  required  to  satisfy 

exactly  the  interlayer  stress  compatibility  conditions  given  by  Eqs.  49,  50a, 

and  50b.  (Note  that  the  condition  on  a is  identically  satisfied  in  this 

z 

case  since  0^=0  everywhere.)  The  conditions  may  be  satisfied  by  appropriate 

choice  of  the  constants  of  integration,  D™  and  D™,  for  each  layer.  Specifically, 
11  12 

and  are  chosen  so  that  Eq.  50a  is  satisfied  (zero  transverse  shear  stress 

at  the  bottom  of  the  laminate).  Then  D™  and  D™  (m=2,3...M)  are  chosen  so  that 

Eq.  49  is  satisfied  at  the  interface  between  layers  m-1  and  m.  This  process 

M M 

is  repeated  until  the  constants  and  have  been  determined.  The  result- 

_ m , m . , 

ing  expressions  for  a and  a are  given  by: 

xz  yz 


<r*7  = < Pa  + * «?  A +«?/&*  * « 7 N 

* k,i  A3  + K!4  A a + KIS  Pt5  + KIS  A*  + Kn  P17  * K!R  P'S> 


V i ' k’fc  A + Ka  f>s  + Ky  / Ks  + Ar  + k& 
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where  the  coefficients  of  the  6*3  are  given  by 
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; the 

z-coordinate  of  the  ith  surface  (see  Fig.  7). 

The  stress  assumption  given  by  Eqs.  83a-83c,  84a,  and  84b  now  satisfies 
the  zero  transverse  shear  stress  requirement  at  the  bottom  surface  of  the 
laminate,  and  the  transverse  shear  stress  compatibility  requirement  at  each 
interlayer  surface  in  the  laminate.  The  only  remaining  requirement  is  that 
the  transverse  shear  stresses  be  zero  on  the  upper  surface  of  the  laminate. 
This  condition  may  be  stated  as 
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(86a) 
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M M 

The  expressions  for  O and  a contain  no  undetermined  constants  of  integra- 

xz  yz 

tion  so  that  Eqs.  86a  and  86b  may  be  viewed  as  constraint  conditions  on  the 
stress  parameters  requiring  elimination  of  two  of  the  8's  in  favor  of  the 
remaining  sixteen  B's.  For  the  present  element,  B^  and  8g  will  be  determined 
in  terms  of  the  remaining  B's.  Thus,  Eqs.  84a  and  84b  are  substituted  into 
Eqs.  86a  and  86b  and  the  resulting  two  equations  are  solved  to  yield,  in 
matrix  form. 
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where  D is  the  determinant  of  the  2 by  2 matrix  coefficient  of  the  vector 

I 8.  80  J and  is  given  by 
4 8 


d - k"  a- 


Kfe  kb 


(88) 


and  where  the  barred  notation  has  been  employed  to  indicate  that  the  coefficient 

of  the  B’s  are  evaluated  at  the  upper  surface  of  tne  laminate  (e.g. 

evaluated  at  z=h  ,).  It  is  important  to  note  that  the  choice  of  which  B's 

to  eliminate  is  not  arbitrary;  they  must  be  chosen  in  such  a way  that  the 

determinant  of  the  coefficient  matrix  of  the  B's  chosen  is  never  zero  (i.e. 

D in  Eq.  88,  or  its  equivalent  if  other  B's  are  chosen  for  elimination). 

An  appropriate  choice  for  the  present  element  is  8^  and  Bg  because  the 

coefficients  and  of  3.  and  8 , respectively,  will  always  be  nonzero. 

The  matrix  multiplication  indicated  in  Eq.  87  can  be  performed  and 

the  resulting  expressions  for  3.  and  8„  can  then,  in  principle,  be  substituted 

4 8 

into  Eqs.  83a,  83b,  83c,  84a,  and  84b  to  yield  a stress  assumption,  in  terms 
of  16  stress  parameters,  which  satisfies  equilibrium  and  which  also  satisfies 
the  interlayer  transverse  shear  stress  compatibility  conditions  and  the 
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requirement  of  zero  transverse  shear  stress  on  the  lower  and  upper  surfaces 

of  the  laminate.  However,  in  practice,  it  is  more  convenient  to  form  the 

laminate  H (Eq.  79a)  and  G (Eq.  79b)  matrices  based  on  the  full  set  of  18  8's 

with  the  matrix  Pm  for  the  mth  layer  being  determined  from  Eqs.  83a-83c,  84a, 

(18)  (18)  (18) 

and  84b  (denote  these  matrices  by  H , G , and  8 ) . A relation  between 

n r j 

the  set  of  18  8's  and  the  final  set  of  16  8's  (denoted  by  8 ) can  be  obtained, 

using  Eq.  87  for  8^  and  8g,  in  the  form 


(89) 


Equation  89  can  be  viewed  as  a transformation  relation  between  sets  of  stress 

parameters;  the  final  laminate  H and  G matrices  (denoted  by  ancj  g^^) 

( 10  ) 

based  on  the  final  set  of  16  stress  parameters,  8 , are  calculated  from: 


H°L)  - HC,&)  T 

~ ~ ~ 18 -lb 


(90a) 
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f/a) 


(90b) 


The  element  B and  k matrices  can  then  be  calculated  from  Eqs.  76  and  77, 
respectively,  using  the  matrices  and 

The  displacement  behavior  of  the  laminate  is  represented  by  transla- 
tional displacements  u,v,  and  w (in  the  x,y,  and  z directions,  respectively) 

of  the  midsurface  of  the  laminate,  and  rotations  9 and  0 of  the  cross- 

x y 

section  normal  lines  about  the  x and  y axes  respectively.  The  rotations 

0 and  9 correspond  to  9w/9y  and  9w/9x  but  are  not  equal  to  these  quantities 
x y 

because  9 and  0 are  assumed  independent  of  w.  Lines  normal  to  the  plate 
x y 

midsurface  prior  to  deformation  need  not  be  normal  to  the  plate  midsurface 

after  deformation.  Thus,  average  laminate  transverse  shear  deformation 

effects  are  included.  The  displacements,  u,v,  and  w along  side  i (between 

nodes  i and  i+1)  are  expressed  in  terms  of  generalized  nodal  degrees  of 

freedom  u,v,w,0  , and  9 (see  Fig.  7)  at  nodes  i and  i+1  by  the  following 
x y 

linear  interpolation 


**•  UL  O-s)  fW(((  S +T*  [ 9y.  (/-s')  + &yt-¥/s] 
tr--  vr-  0-s)  - H [ox.(i-s'>  + 6*.+  (s] 


(91a) 

(91b) 
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W s w,  ( I - S ) 4-  vU 


(91c) 


where  s is  a nondimensional  parameter  taking  on  values  s=0  at  node  i and  s=l 
at  node  i+1.  Note  that  the  inplane  displacements  u and  v vary  linearly  in 
the  z direction,  but  the  transverse  displacement  w is  constant  through  the 
thickness  of  the  laminate.  The  calculation  of  the  matrix  G™  for  the  mth 
layer  is  performed  by  obtaining  the  boundary  displacement  interpolation 
matrix,  L,  from  Eqs.  91,  and  the  boundary  traction  interpolation  matrix, 

Rm,  for  the  mth  layer  from  the  stress  assumption  by 


R"l  -OP"") 

Side  u ^Si'de  l 

(92) 

where 

V is  the  matrix  of  direction  cosines  on  the 

ith  side  of  the  element 

and  is 

obtained  from  the  relations 

’ ^x  + 

(93a) 

= byy  0^  t (S"y 

(93b' 

= ®X2  ^x  +•  ^yz 

(93c) 

where 

-0^  3 COS  OC 

(94a) 

Oy  = Sm  04 

(94b) 

and  where  a is  the  angle  (CCW)  from  the  x-axis  to  the  outward  normal  direction 
for  the  ith  side  of  the  element. 

For  dynamic  analyses,  the  element  mass  matrix  is  required.  Because  of 
the  linear  interpolation  of  displacements  on  the  boundary  of  the  element, 
the  appropriate  choice  for  a displacement  interpolation  in  the  interior  of 
the  present  quadrilateral  element  is  a bilinear  assumption: 
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where  the  IX  terms  are  given  by  Eqs.  72.  The  hybrid-rational  element  mass 

matrix,  m , is  then  obtained  from 
~HR 


MR 


M r 

Z,  J pl  mt#  <=/V 


(96) 
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where  p1  is  the  mass  density  of  the  ith  layer,  V1  is  the  volume  of  the  ith 

n 

layer  of  the  nth  element,  and  where  the  interpolation  matrix,  N,  is  obtained 
from  Eqs.  95.  It  is  important  to  note  again  that  for  the  present  element 
formulation  a reference  surface  z=0  is  chosen  for  the  laminate,  in  contrast 
to  element  ELEMZ  where  a reference  surface  z=0  was  chosen  for  each  layer. 

The  resulting  hybrid-rational  mass  matrix  is  in  general  fully  populated 
and  thus  the  storage  requirements  for  the  assembled  mass  matrix  will  be  the 
same  as  those  for  the  assembled  stiffness  matrix.  An  alternate  approach 
often  used  in  dynamic  analyses  is  to  define  the  mass  properties  of  an  element 
in  terms  of  generalized  nodal  lumped  masses.  In  effect,  the  mass  properties 
of  the  entire  element  are  assumed  to  be  equally  portioned  to  each  of  the 
four  nodes.  The  resulting  lur ped  mass  matrix,  m , is  a diagonal  matrix  (i.e. 

~ Li 

nonzero  contributions  only  on  the  diagonal)  and  is  given  by 
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where 
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and  the  quantities  d^  and  d^  are  given  by 


d,-  4 F P~Um„-0 
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where  A is  the  area  of  the  plate  midsurface,  p is  the  mass  density  of  the 

mth  layer,  h is  the  z coordinate  of  the  mth  interface  (h,=z  coordinate  of 
m 1 

the  bottom  surface  of  the  laminate) , and  M is  the  total  number  of  layers  in 
the  laminate.  The  terms  d^  and  d^  correspond,  respectively,  to  the  transla- 
tional and  rotational  degrees  of  freedom  in  the  element. 

Because  of  the  diagonal  form  of  the  lumped  mass  matrix,  the  assembled 
mass  matrix  can  be  stored  as  a vector  of  length  equal  to  the  total  number 
of  degrees  of  freedom  in  the  assembled  structure,  requiring  significantly 
less  computer  storage  than  the  assembled  mass  matrix  corresponding  to  the 
hybrid-rational  element  mass  matrix.  However,  the  more  important  considera- 
tion in  this  case  is  the  relative  accuracy  of  the  two  approaches  for  predict- 
the  frequencies  and  mode  shapes  of  the  assembled  structure  (as  discussed  in 
Section  7 ) . 

In  subsequent  discussions,  the  present  multilayer  glate  four-node 
quadrilateral  element  will  be  referred  to  as  element  MLP3K(Q).  Corresponding 
to  the  present  development,  computer  subroutines  MLP3K  and  MLP3M  have  been 
written  to  calculate  the  element  stiffness  and  mass  matrices,  respectively, 
and  subroutine  MLP3S  has  been  written  to  calculate  the  stresses  and  strains, 
given  the  calculated  nodal  displacements.  These  modules  can  then  be  utilized 
in  a static  or  dynamic  structural  analysis  package. 

4.3  Modification  to  a Triangular  Element 

Because  of  the  need  to  utilize  refined  mesh  arrangements  in  regions 
of  expected  high  stress  or  strain  gradients  and  more  coarse  mesh  arrangements 
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in  regions  of  the  struccure  where  less  severe  stress  and  strain  gradients 
are  expected,  a transition  element  is  required.  For  present  purposes,  a 
three-node  triangular  element  will  be  used. 

The  requirement  of  displacement  compatibility  on  interelement  boundaries 
dictates  that  the  boundary  displacement  interpolation  used  for  the  triangular 
element  be  identical  to  that  used  for  element  MLP3K (Q) . In  addition,  the 
same  stress  interpolation  used  for  element  MLP3K(Q)  can  be  used  for  the 
triangular  element  (termed  element  MLP3K(T)  in  subsequent  discussions). 

Thus,  the  development  of  element  MLP3K(T)  follows  exactly  the  development 
of  element  MLP3K(Q)  with  the  exceptions  that  the  area  integrals  now  extend 
over  a triangular  planform,  and  that  the  boundary  integrals  now  extend 
over  a three-sided  boundary. 

In  applications  of  the  hybrid-stress  model  to  static  analysis,  the 
accuracy  of  the  element  is  often  dependent  on  how  well  the  stress  and  dis- 
placement assumptions  are  matched.  Because  the  hybrid-stress  model  is  a 
mixed  model,  the  element  becomes  more  stiff  as  the  order  of  the  stress 
assumption  is  increased  (holding  the  displacement  assumption  fixed)  and  the 
element  becomes  more  flexible  as  the  order  of  the  displacement  assumption 
is  increased  (holding  the  stress  assumption  fixed) . Although  no  rigid 
guidelines  are  available  for  determining  the  appropriate  balance  between 
stress  and  displacement  assumptions,  experience  to  date  with  the  quadri- 
lateral multilayer  plate  elements  suggests  that  a reasonable  guideline  is 
that  the  number  of  stress  parameters,  B,  should  be  as  close  to  the  minimum 
number  (corresponding  to  be  3-q  relation  of  Eq.  75)  as  possible.  For 
element  MLP3K(Q)  the  minimum  number  of  B's  is  14  (20  q's  minus  6 rigid-body 
modes)  and  16  B's  are  actually  used  for  that  element;  as  will  be  shown  in 
Section  5,  good  convergence  behavior  is  obtained  for  element  MLP3K(Q).  For 
element  MLP3K (T)  the  minimum  number  of  B's  is  9 (15  q's  minus  6 rigid-body 
modes),  and  16  B's  are  used.  Convergence  studies  given  in  Section  5 will 
show,  as  expected,  that  element  MLP3K (T)  is  significantly  stiffer  than 
element  MLP3K(Q). 

The  present  triangular  element  is  intended  for  use  solely  as  a 
transition  element  and  it  is  expected  that  the  accuracy  of  the  analysis 
will  be  governed  primarily  by  the  accuracy  of  the  quadrilateral  elements 
in  the  mesh.  However,  even  in  the  role  of  transition  element,  the 
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triangular  element  MLP3K(T)  could  have  a degrading  effect  on  the  solution 
by  introducing  bands  of  excessively  stiff  elements  in  the  finite-element 
mesh.  To  alleviate  this  potential  stiffening  effect,  a more  flexible 
triangular  element  is  required.  The  only  alternative  for  obtaining  a 
more  flexible  element  is  to  reduce  the  stress  assumption.  The  displacement 
assumption  cannot  be  increased  because  of  the  displacement  compatibility 
required  between  the  triangular  element  and  element  MLP3K(Q).  Such  an 


element  has  been  developed  in  the  present  study  by  deleting  3, r 3^ r 

4 6 8 9 


3 , and  3. 0 from  the  strain  assumption  of  Eq.  81.  The  resulting  123 

15  18 


stress  assumption  is  obtained  from  Eqs.  83a-83c,  84a,  and  84b  by  deleting 
the  above  B's.  In  order  to  satisfy  the  zero  transverse  shear  stress 
requirements  on  the  upper  surface  of  the  laminate,  3^  and  B-,  are  eliminated 
(these  stress  parameter  numbers  correspond  to  the  original  numbering  used 
in  Eq.  81) . After  reduction,  the  element  contains  10  independent  stress 
parameters,  one  more  than  the  minimum  number  required.  In  subsequent 
discussions,  this  element  will  be  referred  to  as  element  MLTPK. 


4.4  Comparison  of  Alternate  Formulations 

The  fomulation  of  a multilayer  plate  element  (ELEMZ)  applicable  for 
thick-plate  structures  has  been  given  in  Section  3,  and  an  alternate 
formulation  (element  MLP3K (Q) ) for  applications  to  moderately-thick  plate 
structures  has  been  given  in  the  present  section.  The  purpose  of  the 
present  subsection  is  to  compare  the  two  formulations  in  terms  of  applicability, 
computer  core  storage  requirements,  and  computation  time  requirements.  The 
question  of  relative  accuracy  of  the  two  approaches  will  be  addressed  in  the 
next  section. 

In  terms  of  applicability,  it  should  be  expected  that  element  ELEMZ 
will  be  capable  of  adequately  representing  the  severe  cross-sectional  warping 
often  observed  in  thick  laminated  plates,  whereas  element  MLP3K(Q)  will 
represent  this  behavior  only  in  an  average  sense.  Results  presented  by  Mau 
[2]  and  Spilker  [5]  suggest  that  accurate  answers  will  be  obtained  by  using 
element  ELEMZ  for  structures  having  typical  thickness  ratios  (defined  as 
the  ratio  of  the  spanwise  dimension  to  the  thickness  dimension)  as  low  as 
4 (thick  plate),  whereas  the  accuracy  of  element  MLP3K(Q)  will  degenerate 
noticeably  for  aspect  ratios  less  than  8 to  10  (moderately-thick  plate). 
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Most  multilayer  structures  of  interest  will  be  no  more  than  moderately 
thick,  but  will  be  made  up  of  a large  number  of  layers  (i.e.  more  than  20) . 

For  such  structures,  similar  accuracy  should  be  expected  from  both  elements. 
However,  only  element  MLP3K(Q)  will  be  applicable  for  such  problems  in 
practice  because  of  the  excessive  storage  requirements  associated  with 
element  ELEMZ. 

The  clear  superiority  of  MLP3K (Q)  over  ELEMZ  is  found  in  terms  of 
computational  efficiency.  As  discussed  in  Subsection  3.5,  the  storage 
requirements  for  element  ELEMZ  may  be  prohibitively  large  if  the  number  of 
layers  is  large.  This  is  a consequence  of  the  choice  of  a set  of  independent 
stress  parameters  for  each  layer,  and  the  fact  that  the  number  of  nodal 
degrees  of  freedom  depends  on  the  number  of  layers.  For  element  MLP3K(Q) 
the  number  of  stress  parameters  (16)  and  the  number  of  degrees  of  freedom 
(20)  for  an  element  are  fixed  and  thus  the  storage  requirements  for  the 
generation  of  the  MLP3K(Q)  element  stiffness  matrix  and,  on  a global  level, 
for  the  assembled  stiffness  matrix  are  independent  of  the  number  of  layers. 

In  principle,  there  is  no  limitation  on  the  number  of  layers  which  can  be 
accommodated  by  element  MLP3K(Q). 

In  terms  of  computation  time,  element  MLP3K(Q)  will  again  be  superior 

to  element  ELEMZ.  The  generation  of  the  ELEMZ  stiffness  matrix  requires  an 

inversion  of  the  20x20  Hm  matrix  for  each  layer,  whereas  MLP3K (Q)  requires 

only  a single  conversion  of  the  16x16  H matrix  for  the  entire  laminate.  In 

addition,  element  ELEMZ  requires  an  additional  inversion  of  the  matrix  product 
-IT 

AH  A and  the  expression  for  k requires  significantly  more  matrix  multiplica- 
tions than  for  element  MLP3K(Q).  When  element  ELEMZ  is  used,  the  solution 
time  for  the  assembled  matrix  equations  for  the  structure  will  depend  on  both 
the  number  of  elements  and  the  number  of  layers  because  the  total  number  of 
degrees  of  freedom  in  the  assembled  system  depends  on  both.  However,  the 
solution  time  for  element  MLP3K(Q)  will  depend  only  on  the  number  of  elements 
in  the  finite-element  mesh. 

As  discussed  in  Subsection  3.4,  it  is  possible  in  principle  to  develop 
a traction- free-edge  element  based  on  the  formulation  used  for  element  ELEMZ. 
However,  it  appears  that  no  such  traction-f ree-edge  element  can  be  developed 
based  on  the  formulation  used  for  element  MP3K(Q).  It  should  be  recalled 
that  the  stress  assumption  used  for  element  MLP3K(Q)  is  based  on  a fixed  set 
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of  stress  parameters.  However,  the  stresses  0,0,  and  o need  not  be 

x y xy 

continuous  at  interlayer  boundaries.  In  effect,  these  stresses  are  indepen- 
dent from  layer  to  layer  because  the  interpolation  matrix  pm  depends  on 
layer  material  properties  which  vary  from  layer  to  layer.  As  a result,  in 
order  to  satisfy  the  traction-free  condition  in  each  layer,  a reduction  of 
the  total  number  of  B's  would  be  required  for  each  layer  processed.  Clearly, 
this  layer-dependent  reduction  of  B's  is  not  acceptable  for  an  element  in 
which  the  number  of  stress  parameters  is  independent  of  the  number  of  layers. 
Such  an  argument  is  not  limited  to  the  linear  stress  assumption  employed  for 
element  MLP3K(Q);  it  also  applies  to  any  higher  order  assumption. 

Finally,  some  comments  on  the  advantages  of  the  assumed- stress  hybrid 
model  over  the  conventional  assumed-displacement  model  should  be  made.  In 
general,  the  plate  and/or  shell  elements  are  more  easily  formulated  by  the 
hybrid-stress  model  because  compatible  displacement  assumptions  need  be 
made  only  along  the  element  boundaries;  formulations  by  the  assumed-displace- 
ment model  require  the  assumption  of  a displacement  field  in  the  interior 
of  the  element  which  yields  displacement  continuity  along  interelement 
boundaries.  In  many  cases,  assumed-stress  hybrid  elements  yield  improved 
displacement  and  stress  distribution  results  by  comparison  with  assumed- 
displacement  elements  with  similar-order  interpolation.  In  particular,  the 
present  MLP3K(Q)  element  has  been  shown  in  Ref.  5 to  yield  more  accurate 
displacement  results  than  an  assumed-displacement  multilayer  plate  element 
(which  also  includes  transverse  shear  effects)  for  identical  meshes;  however, 
element  MLP3K(Q)  has  five  degrees  of  freedom  per  node  by  comparison  to  the 
seven  degrees  of  freedom  per  node  for  the  assumed-displacement  element. 
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SECTION  5 

ELEMENT  PERFORMANCE  TESTS 


5.1  Cylindrical  Bending  of  a Three-Layer  Infinite  Strip 

The  problem  of  a strip  of  width  £ in  the  y-direction  and  infinite 
length  in  the  x-direction,  of  cross-ply  construction,  total  thickness  h, 
and  loaded  by  a sinusoidal  load  in  the  transverse  direction  has  been  chosen 
as  a first  test  of  the  multilayer  plate  elements.  The  exact  elasticity 
solution  and  the  classical  lamination  theory  solution  have  been  obtained  by 
Pagano  [14].  The  geometry,  material  properties,  coordinate  directions,  and 
loading  function  are  given  in  Fig.  8a.  The  finite  strip  is  assumed  to  be 
simply  supported  along  the  boundaries  y=0  and  y=£. 

Because  the  plate  is  infinitely  long  in  the  x-direction  and  the 
laminate  constuction  is  balanced,  the  finite-element  analysis  may  be  per- 
formed by  analyzing  a strip  of  finite  width,  a,  in  the  x-direction,  and  of 
length  £/2  because  of  symmetry  in  the  y-direction  (see  Fig.  8b).  Ten  equally 
sized  square  elements,  with  side  length  a=£/20,  are  employed  in  the  finite- 
element  analysis.  Simply-supported  boundary  conditions  are  imposed  on  the 
side  y=0,  and  symmetry  conditions  are  imposed  on  the  other  three  sides. 

Elements  ELEMZ  and  MLP3K(Q)  are  employed  in  the  present  analysis  to  obtain 
a comparison  of  their  performances. 

The  following  quantities  of  interest  have  been  chosen  to  obtain  a 
comparison  of  the  results  obtained  by  using  each  of  the  above  elements:  (1) 

the  normal  stress  a at  y=£/2;  (2)  the  transverse  shear  stress  x at  the 

y yz 

boundary  y=0;  (3)  the  inplane  displacement  v in  the  y-direction  at  the  boundary 

y=0;  and  (4)  the  midplane  transverse  displacement  w in  the  z-direction  at  the 
center  of  the  plate  y=£/2.  For  convenience,  the  results  will  be  presented  in 
terms  of  normalized  quantities  (barred)  which  are  given  by 


H I 


where  is  the  amplitude  of  the  sinuisoidal  loading  and  E ^ is  the  elastic 


modulus  perpendicular  to  the  fiber  direction.  It  should  be  noted  that  the 
quantity  S is  a measure  of  the  thickness  ratio  of  the  plate  (S=4  is 


considered  to  be  a thick  plate,  S=50  is  considered  to  be  a thin  plate) . 


Results  are  presented  for  a three-layer  (90°/0°/90°)  laminate  (N=3) 


and  a seven- layer  (90O/0O/90O/0°/90O/0O/90O)  laminate  (N=7) . Two  subcases 


are  then  considered  for  each  of  the  above  laminates:  one  in  which  the  plate 
is  thick  (S=4)  and  one  in  which  the  plate  is  moderately  thick  (S=10) . For 
each  case,  and  in  the  corresponding  tables  and  figures,  results  are  presented 
which  were  obtained  by  an  exact  solution  [14],  classical  lamination  theory 
[14],  and  by  using  elements  ELEMZ  and  MLP3K (Q ) . Table  1 qives  the  results 
obtained  for  the  transverse  displacement,  w,  at  y->' /2.  The  results  obtained 
by  using  element  ELEMZ  are  in  excellent  agreemen  with  the  exact  solution 
for  all  cases  (0.8%  error  for  S=4,  N=3,  0.003%  error  for  S=10,  N=3,  0.63% 
error  for  S=4,  N=7,  and  0.15%  error  for  S=10,  N=7).  The  results  obtained 
by  using  element  MLP3K (Q)  for  the  three-layer  cases  are  too  flexible  for 
the  case  S=4  (9.51%  error),  but  are  in  excellent  agreement  with  the  exact 
solution  for  the  case  S=10  (0.96%  error) . The  results  obtained  by  using 
element  MLP3K(Q)  for  the  seven  layer  laminate  are  in  good  agreement  with 
the  exact  solution  for  both  S=4  (1.81%  error)  and  S=10  (0.08%  error),  which 
suggests  that  for  increasing  number  of  layers  the  average  cross-sectional 
rotation  behavior  incorporated  in  element  MLP3K(Q)  more  closely  models  the 
exact  behavior.  Note  that  the  lamination  theory  solution  is  in  poor  agree- 
ment with  the  exact  solution  even  for  moderately-thick  laminates  (S=10) 
which  suggests  that  transverse  shear  effects  (neglected  in  lamination  theory) 
play  a significant  role  in  the  calculation  of  transverse  deflection. 


The  results  obtained  for  the  stresses  O and  x for  the  three-layer 

y yz 


cases  are  plotted  in  Figs.  9a  and  9b  as  a function  of  the  nondimensionalized 
thickness  parameter  z.  Note  that  the  results  are  plotted  only  from  z=0  to 


0.5  because  the  quantities  a and  T are  odd  and  even  functions  of  z, 

y yz  _ 

respectively.  As  is  shown  in  Fig.  9a,  the  results  obtained  for  by  using 


element  ELEMZ  are  in  excellent  agreement  with  the  exact  solution,  whereas 
the  results  obtained  by  using  element  MLP3K(Q)  are  in  essential  agreement 
with  the  classical  lamination  theory  solution.  Note  that  the  lamination 
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theory  solution  gives  poor  agreement  with  the  exact  solution  for  S=4,  but 
much  better  agreement  for  S=10. 

The  results  for  the  transverse  shear  stress,  T for  the  three-layer 

yz 

laminate  (Fig.  9b)  obtained  by  using  element  ELEMZ  are  in  reasonable  agree- 
ment with  the  exact  solution  for  both  S=4  and  S=10,  although  the  tendency 
appears  to  be  to  underestimate  slightly  the  interlaminar  shear  stress  value 
at  the  interface  between  the  top  and  middle  layers,  and  to  overestimate 
slightly  the  value  at  the  center  (z=0)  of  the  laminate.  The  results  obtained 
by  using  element  MLP3K(Q)  are  in  general  agreement  with  the  lamination  theory 
solution  in  the  top  layer,  but  tend  more  toward  the  exact  solution  in  the 
center  layer.  For  the  case  S=10,  the  MLP3K (Q)  results  are  actually  slightly 
more  accurate  than  the  ELEMZ  results. 

The  inplane  displacement  v,  given  in  Fig.  9c,  demonstrates  the  degree 
of  cross-sectional  warping  present  in  the  three-layer  laminate.  For  the 
case  S=4,  severe  warping  is  present,  but  for  the  case  S=10,  only  moderate 
warping  is  present.  For  both  cases,  the  results  obtained  by  using  element 
ELEMZ  are  in  excellent  agreement  with  the  exact  solution,  whereas  the  results 
obtained  by  using  element  MLP3K(Q)  are  in  essential  agreement  with  the  lamina- 
tion theory  solution.  The  lamination  theory  solution  gives  a good  approxima- 
tion to  the  exact  solution  for  S=10,  but  a poor  approximation  to  the  exact 
solution  for  S=4. 

The  results  obtained  for  the  stresses  0 and  T and  the  inplane  dis- 
_ y yz 

placement,  v,  for  the  seven- layer  laminate  are  shown  in  Figs.  10a  through 
10c.  The  observations  made  about  the  element  behavior  for  the  three-layer 
laminate  also  hold  for  the  seven-layer  laminate.  The  results  obtained  by 
using  element  ELEMZ  are  found  to  be  in  good  agreement  with  the  exact  solution 
for  all  cases,  and  are  clearly  superior  to  the  results  obtained  by  using 
element  MLP3K (Q)  for  the  cases  corresponding  to  a thick  laminate  (S=4) . For 
the  moderately  thick  cases  (S=10) , the  accuracy  of  the  MLP3K(Q)  results  is 
comparable  to  or  slightly  better  than  the  ELEMZ  results. 

For  the  present  ten-element  (3-layer)  finite  element  solutions,  approxi- 
mately 215,000  BYTES  of  computer  core  storage  (in  double  precision)  were 
required  for  the  ELEMZ  solutions,  whereas  only  106,000  BYTES  of  computer  core 
storage  were  required  for  the  MLP3K(Q)  solution.  For  the  seven- layer  laminate 
cases,  using  the  same  10  element  mesh,  430,000  BYTES  of  storage  were  required 


for  the  ELEMZ  solutions,  whereas  the  MLP3K(Q)  solution  still  required  only 
106,000  BYTES  of  storage.  When  the  4 cases  (3  layers,  7 layers,  S=4,  and 
S=10)  were  run  together,  the  ELEMZ  solution  required  approximately  0.461 
CPU  minute,  and  the  MLP3K(Q)  solution  required  approximately  0.040  CPU 
minute.  Thus,  significant  savings  in  computer  core  storage  and  in  computa- 
tion time  were  found  when  using  element  MLP3K(Q). 

A comparison  of  the  results  obtained  by  using  elements  ELEMZ  and 
MLP3K(Q)  shows  element  ELEMZ  to  be  (as  expected)  the  more  accurate  element, 
particularly  for  the  case  S=4  (thick  laminate).  However,  for  the  case  S=10, 
the  superiority  of  element  ELEMZ  over  element  MLP3K(Q)  becomes  much  less 
significant,  and  the  MLP3K(Q)  results  are  well  within  the  accuracy  limits 
of  engineering  analysis.  In  view  of  the  large  storage  requirements  associated 
with  ELEMZ  as  the  number  of  layers  is  increased,  it  may  be  expected  that  the 
use  of  ELEMZ  must  be  restricted  to  cases  where  the  number  of  layers  is  small 
and  the  laminate  is  thick  (e.g.  the  case  S=4) . For  cases  where  a representa- 
tive thickness  aspect  ratio  such  as  S is  greater  than  10,  element  MLP3K(Q) 
can  be  used  with  confidence  regardless  of  the  number  of  layers.  Based  on 
the  results  of  the  present  example,  it  may  be  expected  that  if  applied  to  thick 
laminated  plate  problems,  element  MLP3K(Q)  will  yield  transverse  normal  dis- 
placement results  which  are  significantly  more  accurate  than  lamination  theory, 
all  other  results  tending  toward  lamination  theory.  However,  it  should  be 
recalled  that  element  MLP3K(Q)  is  not,  strictly  speaking,  based  on  lamination 
theory.  Transverse  shear  stresses,  strains,  and  strain  energy  are  included 
automatically  in  element  MLP3K(Q),  but  transverse  shear  strain  energy  is  not 
included  in  lamination  theory.  Finally,  for  thick  laminates  of  few  layers 
(e.g.  S=4,  N=3) , it  appears  that  the  accuracy  of  the  stress  distribution 
through  the  thickness  is  most  strongly  influenced  by  the  cross-sectional 
warping  behavior  of  the  laminate. 

5.2  A Square  Two- Layer  P 1 ate  Under  Uniform  Transverse  Loading 

The  example  problem  to  be  considered  here  is  a square  plate  of  side 

length,  a,  of  two- layer  angle-ply  construction  (+0  fiber  orientations  with 

respect  to  the  global  x axis) , loaded  by  a uniform  pressure  load  of  magnitude 

q . The  plate  has  total  thickness  h,  and  each  layer  is  of  equal  thickness, 
o 

h/2.  The  plate  is  simply  supported  on  all  four  edges  such  that  displacements 
transverse  and  normal  to  the  edge  are  not  permitted,  but  motion  in  a direction 
parallel  to  the  boundaries  is  permitted.  Each  layer  is  composed  of  the  same 
material;  the  material  properties,  dimensions,  and  geometry  are  given  in  Fig.  11. 
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A solution  for  the  present  problem  based  on  classical  lamination  theory 

has  been  presented  by  Whitney  [15]  using  a Fourier  series  approach.  The 

equations  required  to  calculate  the  Fourier  coefficients,  displacements, 

stress  resultants,  and  moment  resultants  are  given  in  Ref.  15.  However,  the 

expressions  given  for  the  moment  resultants,  M and  M , appear  to  be  in  error. 

x y 

For  convenience,  the  correct  expressions  for  M and  M will  be  derived  here. 

x y 

The  inplane  displacements,  u and  v , at  the  laminate  midsurface,  and 
transverse  displacement,  w,  are  expressed  as  a Fourier  series  in  the  form 
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where  a and  b are  the  dimensions  of  the  plate  in  the  x and  y directions, 

respectively,  and  the  Fourier  coefficients,  E , F , and  G are  obtained 
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by  expressions  given  in  Ref.  15.  From  the  laminate  constitutive  relations. 
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where  the  A^ . , B _ , and  D _ terms  are  given  for  the  present  two-layer 
laminate  by 
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and  where  C. . are  the  reduced  (plane-stress)  stiffness  material  properties 
* 1 

of  the  mth  layer  (Eqs.  42)  and  h and  h , are  the  z-coordinates  of  the 

m m-1 

upper  and  lower  surface  of  the  mth  layer.  For  +0  angle-plies,  it  can  be 

shown  that  the  terms  B.,,B,_,B  and  D„ , are  zero,  so  that  the  expressions 

11  12  22  16  26 

for  the  moment  resultants,  M and  M , are  obtained  from  Eq.  102  as 

x y 
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The  necessary  inplane  strain  and  curvature  terms  can  be  obtained  by 
differentiation  of  Eqs.  101a  through  101c: 
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Substituting  Eqs.  105a  through  105c  into  Eqs.  104a  and  104b  yields  the  final 

expression  for  M and  M : 
x y 
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(106c) 


where  R is  the  spanwise  aspect  ratio  and  is  given  by 


o. 


The  exact  solutions  for  the  present  example  problem  have  been  numerically 
generated  by  using  Whitney's  expression  (except  for  and  where  Eqs. 

106a  and  106b  have  been  substituted)  and  retaining  400  terms*  in  the  Fourier 
series . 

Because  the  results  presented  in  the  previous  subsection  suggest  that 
element  MLP3K(Q)  may  be  the  more  widely  useful  element,  the  present  example 
problem  has  been  chosen  to  obtain  additional  performance  information  for 
element  MLP3K(Q).  Because  of  the  bending/extensional  coupling  in  the  present 
unbalanced  laminate,  symmetry  conditions  cannot  be  invoked,  and  the  entire 
plate  must  be  modeled  in  the  finite-element  analysis.  A K-by-K  mesh  of 
square  elements  is  employed,  where  K is  the  number  of  elements  in  each 
direction. 

The  first  case  to  be  considered  is  a plate  of  side  length  a=10",  thick- 
ness h=0.2",  fiber  orientation  angles  9=+45°,  and  a uniform  load  qo=100  psi. 

The  quantities  of  interest  are:  (1)  the  transverse  deflection  w at  the  center 

of  the  plate;  (2)  the  in-plane  x-direction  displacement  u°(z=0)  at  x=a/2,  y=0; 
and  (3)  the  moment  at  the  center  of  the  plate.  The  percent  errors  (compared 
to  Whitney's  Fourier  series  approach)  in  the  finite-element  results  (using 
element  MLP3K(Q))  for  these  quantities  are  plotted  in  Fig.  12a  versus  the 
total  number  of  degrees  of  freedom  in  the  finite-element  solution.  The 
points  correspond,  respectively,  to  values  of  K=4,6,8,10,  and  12  (K-by-K 
mesh) . The  results  for  u°  converge  very  rapidly  to  the  exact  solution.  The 
convergence  of  the  results  for  w is  much  slower,  but  for  the  12-by-12  mesh, 
the  w solution  is  in  error  by  only  1.5%.  The  results  for  are  converging 
rapidly,  but  have  not  yet  reached  the  exact  solution  for  the  12-by-12  mesh 
(0.7%  error).  Because  the  stress  (and,  thus,  moment)  solution  is  dependent 
upon  both  in-plane  and  out-of-plane  displacements  for  the  present  unbalanced 
laminate,  this  solution  should  not  converge  to  the  exact  solution  until  both 
u°  and  w have  converged.  Note  that  for  all  cases,  the  moment  is  calculated 
from  the  finite-element  solution  by  integrating  the  average  stress  distribution 


★ 

Significantly  more  terms  than  required  to  obtain  a converged  solution. 
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(obtained  by  nodal  averaging  at  the  center  of  the  plate) . Note  also  that 
for  the  present  example  (9=+45°) , u°  at  x=a/2,  y=0  is  the  same  as  v°  (in- 
plane displacement  in  the  y-direction)  at  x=0,  y=a/2,  and  M =M  at  the 

x y 


center  of  the  plate. 

To  assess  the  effects  of  the  fiber  orientation  angle,  0,  on  the  finite- 

element  solution,  a second  case  has  been  considered.  In  this  case  the  plate 

dimensions  and  loading  are  identical  to  the  first  case  (a=10",  h=0.2”,  and 

qQ=100  psi) , and  a 10-by-10  mesh  has  been  used  to  model  the  entire  plate.  A 

series  of  results  was  obtained  for  the  fiber  orientation  angles,  +0=5°,  15°, 
o o o 

25  , 35  , and  45  . The  quantities  of  interest  for  the  present  case  are  those 
defined  for  the  first  case  and,  in  addition,  the  in-plane  displacement  v°  in 
the  y-direction  at  x=0,  y=a/2,  and  the  moment  at  the  center  of  the  plate. 
The  percent  errors  in  the  finite-element  solution  for  these  quantities  versus 
the  fiber  orientation  angle,  +0,  are  plotted  in  Fig.  12b.  The  exact  solution 
finite-element  solution,  and  error  in  the  finite-element  solution  for  each 
of  these  values  of  0 are  also  tabulated  in  Table  2.  Figure  12b  shows  that 
the  absolute  value  of  the  error  for  these  quantities  is  less  than  3 percent 
for  all  values  of  0.  The  mild  fluctuations  in  the  error  as  0 is  varied  may 
be  attributed  to  the  increase  (or  decrease)  in  the  relative  significance  of 
a particular  quantity  (see  Table  2).  For  example,  a comparison  of  the  exact 
solutions  for  u°  and  v°  shows  that  u°  is  small  compared  with  v°  for  +0=5° 
and  15°,  and  the  finite-element  solution  for  u°  is  less  accurate  until  the 
values  of  u°  and  v°  are  quite  close.  A similar  observation  can  be  made  for 
and  M^.  This  kind  of  behavior  is  typical  of  most  numerical  procedures 
requiring  the  solution  >f  a large  number  of  simultaneous  equations;  the 
dominant  terms  are  predicted  more  accurately.  However,  it  should  be  noted 
again  that  the  fluctuations  in  the  error  for  the  quantities  of  interest  as 
0 is  varied  are  quite  small  and  all  quantities  of  interest  are  predicted  to 
within  3 percent  error  regardless  of  the  value  of  0. 

The  final  performance  test  for  element  MLP3K(Q)  is  the  effect  of  element 
spanwise  aspect  ratio.  For  this  case,  a rectangular  plate  of  two-layer  +45° 
angle-ply  construction  is  considered  with  material  properties  and  loading 
identical  to  case  1.  The  plate  dimension  in  the  y direction,  b,  is  fixed 
at  10"  and  the  plate  dimension  in  the  x direction,  a,  is  increased  to  give 
values  of  the  spanwise  aspect  ratio,  R (Eq.  106c),  of  1, 2, 3, 4 , 5 , 6, 8 , and  10. 
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A 10x10  mesh  is  used  for  all  values  of  R to  model  the  entire  structure,  so 

that  in  all  cases  R is  also  the  spanwise  aspect  ratio  of  each  element. 

Shown  in  Fig.  12c  are  the  percent  errors  in  the  quantities  u°  at  x=a/2, 

y=0,  v°  at  x=0,  y=b/2,  and  w,  M , and  M at  the  center  of  the  plate  (x=a/2, 

x y 

y=b/2)  as  functions  of  R.  As  shown,  the  effect  of  increased  aspect  ratio 

is  to  stiffen  the  predicted  response.  A very  mild  stiffening  effect  is 

observed  for  the  quantities  u°,  w,  M , and  M , with  very  little  degeneration 

x y 

of  the  solution  even  for  the  extreme  case  of  R=10.  However,  a severe 
stiffening  effect  is  observed  for  the  quantity  v°.  The  increasing  error 
in  the  quantity  v°  may  be  viewed  largely  as  a modeling  error  rather  than 
ill-conditioning  of  the  element  as  a function  of  aspect  ratio.  Fixing  the 
number  of  elements  in  the  x direction  and  then  increasing  the  plate  dimension 
in  the  x direction  is  analogous  to  fixing  the  dimension  and  then  decreasing 
the  number  of  elements,  and  the  effects  of  either  approach  are  the  same. 

In  effect,  fewer  elements  are  available  to  model  the  distribution  of  v° 
along  the  x direction  of  the  plate.  In  contrast,  the  dimension  in  the 
y direction  and  number  of  elements  in  the  y direction  are  fixed  in  -the 
present  example  and  the  percent  error  in  the  quantity  u°  at  x=a/2,  y=0 
increases  only  slightly  as  R increases  (i.e.  the  accuracy  in  the  modeling 
of  the  distribution  of  u°  along  the  y direction  is  nearly  independent  of  R) . 
Difficulties  associated  with  the  use  of  large  element  aspect  ratio  are 
usually  manifested  by  severe  ill-conditioning  of  the  assembled  matrix  equa- 
tions resulting  in  large  numerical  errors  in  all  quantities  of  interest. 

No  such  difficulties  have  been  observed  for  the  present  cest  problem  which 
utilizes  rectangular  elements  only.  However,  in  an  application  example 
described  in  Subsection  8.2,  the  H matrix  for  a quadrilateral  element  of 
spect  ratio,  R=7,  could  not  be  inverted.  Because  of  this,  the  use  of 
element  aspect  ratios  greater  than  5 is  not  recommended. 

Limited  results  have  been  obtained  by  using  the  triangular  shap<=>8 
elements  MLP3K(T)  and  MLTPK.  The  example  chosen  is  identical  to  the  first 
case  presented  in  this  subsection  (i.e.  2 layers,  0=+45°,  square  plate). 

The  entire  plate  is  modelled  with  a K-by-K  uniform  mesh  where  each  square 
< ontains  two  triangular  elements.  The  results  of  a convergence  study  using 
lement  ML  3K (T)  are  shown  in  Fig.  12d  for  K=4,6,8,  and  10.  The  finite 
element  solution  for  the  quantities  of  interest  (v°  at  x=0,  y=a/2,  w at  the 


center  of  the  plate,  and  at  the  center  of  the  plate)  is  in  error  by  more 
than  35%  even  for  the  10x10  mesh  (200  triangular  elements) . This  element 
is  clearly  unacceptable  for  use  in  modeling  the  entire  structure,  and  is 
also  felt  to  be  too  stiff  for  use  as  a transition  (mesh  expander)  element. 

It  should  be  recalled  that  element  MLP3K (T)  is  based  on  the  same  stress 
assumption  used  in  element  MLP3K(Q).  The  results  of  a convergence  study 
using  element  MLTPK  (based  on  a reduction  of  the  stress  assumption  used  in 
element  MLP3K(Q))  are  shown  in  Fig.  12e.  Although  the  convergence  behavior 
is  not  monotonic,  it  is  substantially  better  than  that  observed  for  element 
MLP3K (T) . The  use  of  element  MLTPK  to  model  the  entire  structure  is  not 
recommended,  but  this  element  is  suitable  for  use  as  a transition  element. 
Consequently,  element  MLTPK  should  be  used  for  all  multilayer  plate/shell 
problems  requiring  a transition  element. 

The  example  problems  in  the  present  subsection  have  been  chosen  to 
verify  the  accuracy  and  range  of  applicability  of  element  MLP3K(Q).  The 
application  of  element  MLP3K(Q)  to  the  static  analysis  of  a finite  rectangular 
plate  with  a circular  cutout  subjected  to  four-point  bending  is  discussed  in 
Subsection  8.2.  Additional  performance  tests  of  element  MLP3K(Q)  for  dynamic 
analyses  are  given  in  Subsection  8.5. 


SECTION  6 


INTEGRALLY-STIFFENED  PLATES  AND  SHELLS 


6.1  Introduction 

Stiffeners  are  of  interest  in  plates  and  shells  since  they  are  commonly 
used  to  reduce  stress  concentrations  in  panels  around  cutouts  or  to  increase 
the  buckling  strength  of  webs.  There  are  two  ways  to  model  stiffened  plates. 

The  first  approach  is  the  nonintegrally  stiffened  method  which  divides 
the  stiffeners  and  plate  into  separate  elements,  using  beam  elements  to 
model  the  stiffeners.  Both  elements  are  easy  to  formulate,  but  sometimes 
displacement  incompatibility  exists  between  the  stiffener  and  plate  inter- 
faces. This  occurs  in  the  present  study,  which  uses  the  quadrilateral 
element  developed  in  Section  4 and  beam  elements  based  on  the  displacement 
model.  The  boundary  transverse  displacement  for  the  quadrilateral  plate 
element  is  linear  and  this  is  incompatible  with  the  cubic  interpolation  of 
the  beam  stiffener.  This  approach  is  also  inefficient  when  the  stiffeners 
are  closely  spaced,  since  nodal  lines  must  fall  along  stiffeners  and  hence 
the  number  of  plate  elements  increases. 

The  second  approach  is  the  integrally-stiffened  method  which  avoids 
some  of  the  difficulties  mentioned  above  by  formulating  the  plate-stiffener 
combination  as  a single  element.  Nodal  lines  can  then  be  spaced  conveniently 
without  regard  to  stiffener  locations  and  displacement  compatibility  is 
satisfied.  A formulation  for  an  integrally-stiffened  element  based  on  the 
displacement  model  is  developed  in  Ref.  16.  An  equivalent  element  based 
on  the  hybrid  model  is  developed  in  Subsection  6.2.  In  the  latter  formula- 
tion it  is  difficult  to  satisfy  stress  compatibility  at  the  stiffener- 
plate  interface.  Also,  the  stress-free  conditions  are  not  satisfied  at  the 
stiffener  sides  and  at  one  free  surface  of  the  plate. 

A test  problem  of  an  isotropic  simply-supported  stiffened  plate  under 
uniform  transverse  loading  is  analyzed  by  the  above  two  methods.  These 
results  are  compared  in  Subsection  6.3  with  a Ritz  analysis  of  the  plate. 

Kirk  [17]  has  shown  that  the  Ritz  method  is  extremely  accurate,  and,  thus, 
it  provides  an  independent  standard  for  assessment  of  the  finite-element 
methods . 


6.2  Integrally-Stiffened  Quadrilateral  Element 

The  integrally-stiffened  multilayer  plate  element  is  based  on  the 
hybrid  model  and  follows  the  same  basic  formulation  as  the  quadrilateral 
element  MLP3K(Q)  discussed  in  Subsection  4.2.  The  stiffener  is  treated  as 
a finite  layer  attached  to  the  plate  (Fig.  14)  and  the  following  approxima- 
tions are  made: 

1.  The  hybrid  element  is  based  on  the  Principle  of  Modified  Comple- 
mentary Energy  which  requires  that  the  stresses  be  in  equilibrium 
in  the  element.  Thus,  stress  compatibility  must  be  satisfied  at 
layer  interfaces.  This  is  done  in  the  present  MLP3K (Q)  element 
with  stress  assumptions  that  include  transverse  shear  stresses  which 
vary  only  through  the  plate  thickness.  However,  at  the  plate-stiffener 
interface,  which  represents  only  a portion  of  the  plate  area,  the 

fact  that  the  transverse  shear  stresses  are  constant  over  the  plate 
surface  implies  that  only  shear  force  rather  than  stress  compatibility 
is  satisfied,  i.e.,  the  equilibrium  condition  is  met  in  an  integral 
sense  rather  than  exactly.  The  above  approximation  is  chosen  to 
be  consistent  with  the  present  MLP3K(Q)  stress  assumption.  It  is 
not  the  only  possible  choice,  but  other  assumptions  could  not  be 
investigated  within  the  scope  of  the  present  study.  The  stresses 
assumed  also  violate  certain  stress-free  boundary  conditions.  The 
stress-free  condition  on  the  surface  of  the  plate  outside  of  the 
stiffener  is  violated.  On  the  stiffener  itself  the  stress-free 
conditions  on  the  lateral  sides  are  violated.  (However,  note  that 
these  stress-free  conditions  apply  to  surfaces  which  are  not  a 
part  of  the  inter- layer  boundary,  and  therefore,  they  are  not 
required  to  be  satisfied  by  the  modified  complementary  energy 
principle. ) 

2.  The  entire  stiffener  property  is  lumped  on  a plane  perpendicular 
to  the  plate.  Hence,  no  variation  of  stress  across  the  width  of 
the  stiffener  is  considered. 

3.  It  is  assumed  that  the  neutral  surface  of  the  plate  remains 
unaffected  by  the  stiffener. 


4.  Stiffeners  are  modelled  on  the  lower  surface  of  the  plate  only 

(Fig.  14).  This  is  an  arbitrary  restriction;  equivalent  formula- 
tions are  possible  with  stiffeners  on  the  upper  surface  or  on  both 
surfaces. 

These  approximations  limit  the  stiffened  plate  element.  Since  the 
stiffener  is  lumped,  the  element  is  restricted  to  stiffeners  which  are 
narrow  relative  to  the  dimensions  of  the  plate  element.  The  stiffeners 
are  also  restricted  to  small  moments  of  inertia,  such  that  the  neutral 
surface  of  the  plate  is  not  significantly  altered. 

The  strain  assumption  used  in  the  quadrilateral  element  (Eq.  81)  is 
also  assumed  in  the  integrally-stiffened  plate  element.  However,  the 
resulting  equations  for  transverse  shear  stresses  and  the  H and  G matrices 
(Subsection  4.2)  change  slightly  as  a result  of  the  above  approximations. 

Since  the  stiffeners  are  treated  as  layers,  Eqs.  79  become 

ti  = Hp  -f~  H p -h  + Hp  i-  )-)^  -h  ~h  H5  (107a) 
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where  the  subscript  'p'  refers  to  the  plate  and  's’  to  the  stiffeners  and 
where 


& 


(108b) 


The  boundary  surface  9V1  (Fig.  14)  includes  only  the  interelement  boundaries 

s 

of  the  ith  stiffener. 

The  transverse  shear  stresses  exactly  satisfy  the  stress-free  condition 
at  the  bottom  surface  of  the  stiffener.  However,  the  constraints  on  the 
shear  stresses  change  slightly  from  Subsection  4.2  due  to  the  approximation 
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of  shear  force  compatibility.  The  shear  stress  in  the  ith  stiffener  is 
given  by 
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Except  for  the  foregoing  changes  in  the  transverse  shear  stresses  and  the 
H and  G matrices,  the  remaining  formulation  for  the  quadrilateral  element 
stiffness  matrix  follows  that  of  Subsection  4.2. 


6.3  Example  Problem  and  Results 


The  example  problem  is  a simply  supported  square  isotropic  plate 
stiffened  symmetrically  by  four  stiffeners  (Fig.  13).  The  plate  is  loaded 
by  a uniform  transverse  pressure  p^.  The  problem  is  analyzed  by  using 
three  methods. 


Since  the  plate  has  two  axes  of  symmetry  only  a quarter  of  the  plate 
will  be  considered  for  finite-element  analysis.  In  the  first  analysis, 
the  plate  is  modelled  by  element  MLP3K(Q)  developed  in  Section  4 and  the 
stiffeners  are  modelled  separately  by  a beam  element  based  on  the  displace- 
ment model  with  shear  effects  included.  The  displacement  interpolation  in 
the  beam  is  cubic  for  bending  and  linear  for  torsion  and  extension.  It 
has  two  nodes  and  six  degrees  of  freedom  per  node  (three  translations  and 
three  rotations).  Since  the  boundary  displacements  of  the  plate  element 
are  linearly  interpolated  between  nodes,  the  transverse  displacements  will 
be  incompatible.  Four  different  meshes  are  used:  2x2,  4x4,  6x6,  and  8x8. 
These  same  meshes  are  again  modelled  in  the  second  analysis  which  uses  the 
integrally  stiffened  element  developed  in  Subsection  6.2. 

The  third  analysis  applies  the  Ritz  method  by  modelling  the  deflection 
behavior  of  the  plate  as 


here  'a'  is  the  length  of  the  plate.  The  assumed  displacement  field 
satisfies  the  prescribed  force  and  displacement  boundary  conditions  for  a 
simply-supported  square  plate. 


The  total  potential  energy  of  the  structure  can  be  evaluated  from  the 
assumed  deflection  in  terms  of  the  unknown  coefficients  A^,  and  A . 

The  minimization  of  the  potential  energy  with  respect  to  these  coefficients 
then  gives  a set  of  simultaneous  linear  equations  in  these  unknowns  which 
can  be  solved  to  evaluate  them. 

The  total  potential  energy  tt^  of  the  structure  can  be  written  as 


TT  = U-W 
P 


where 

U = strain  energy  of  the  structure 
W = work  done  by  applied  loads 

The  strain  energy  can  be  evaluated  for  the  plate  and  stiffener  and  the 
total  strain  energy  is  written  as 

U = U + U 
s p 


(114a) 


(114b) 
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where 


(115b) 


Under  a uniform  transverse  pressure  loading  the  work  done  is  given  by 


(115c) 
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since  the  plate  is  thin,  transverse  shear  deformation  effects  have  been 
neglected. 
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The  transverse  center  deflection  of  the  plate  evaluated  by  the  three 

methods  has  been  tabulated  in  Table  3 for  comparison.  It  is  apparent  that 

the  nonintegral ly  stiffened  method  is  acceptably  accurate.  The  center 

deflections  computed  by  this  method  are  slightly  greater  than  those  computed 

by  the  Ritz  method  since  the  finite  element  formulation  includes  transverse 

shear  deformation  effects.  However,  for  the  integrally-stiffened  element, 

the  error  is  about  14  percent  with  t =0.1  in.  and  w = 0.03  in.  (Fig.  13). 

s s 

For  the  2x2  mesh  the  area  ratio  AR  = 0.01  and  with  an  8x8  mesh,  AR  = 0.04 
(the  area  ratio  is  the  same  for  both  stiffeners) . Thus,  the  nonintegrally 
stiffened  method  has  provided  an  accurate  analysis,  but  the  integrally- 
stiffened  element  is  poor  and  hence  this  particular  formulation  should  be 
rejected. 

In  any  case,  it  must  be  recognized  that  the  use  o„  beam-like  behavior 

in  the  stiffener  model  restricts  the  analysis  to  the  regime  of  small  stiffeners, 

i.e.  with  area  ratio  AR<0.1.  Stiffeners  which  exceed  this  size  may  begin  to 

have  a plate-like  influence  upon  the  structure  in  either  of  two  ways.  If 

w is  large,  then  the  structure  should  really  be  treated  as  an  assembly  of 
s 

plates  having  different  thicknesses.  On  the  other  hand,  if  tg  is  large, 
then  a three-dimensional  model  is  required,  with  the  stiffeners  represented 
by  plate  elements  perpendicular  to  the  primary  structure. 
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SECTION  7 
DYNAMIC  ANALYSIS 


7.1  Introduction 

New  modules  for  dynamic  analysis  were  added  to  the  existing  static- 
analysis  FEABL-2  software  [6]  during  the  present  investigation.  The  new 
software  consists  of  a collection  of  subroutines  for  solution  of  the 
undamped  free- vibration  eigenvalue  problem: 

( K - ' 0 (117) 

and  one  subroutine  which  combines  the  eigen-solutions  with  load  time 
histories  for  transient  resporse  calculations.  The  eigen-analysis  procedure 
is  based  upon  the  same  approach  used  in  an  earlier  computer  code  [2,3],  but 
it  has  been  re-programmed  to  achieve  modularization  and  hardware- independence. 
The  transient  response  procedure  has  been  programmed  to  permit  analyses  of 
damped  structures  (a  feature  not  included  in  the  earlier  code) , as  well  as 
undamped  structures. 

During  development  of  the  above  software,  it  was  determined  that  some 
modifications  were  required  in  the  existing  FEABL-2  code  in  order  to  achieve 
inter-subroutine  compatibility.  These  modifications  were  made  and  were 
exercised  as  part  of  the  software  verification  program  (Subsection  7.5.3). 

The  resulting  modified  software  has  been  designated  FEABL-5*  and  will  be 
documented  separately  in  a self-contained  user's  guide.  Only  the  numerical 
analysis  methodology,  general  software  organization,  and  verification  tests 
are  reviewed  in  this  report. 

7.2  Eigenvalue  Solution  Method 

The  f ree-vibration  eigenvalue  problem  is  posed  in  the  following  form 
for  numerical  solution: 

Ku--[iy  x (ns) 

where 

y * ( fei}  ■"  tej ) <n9) 

* 

The  capabilities  of  the  FEABL-5  software  include  both  static  and  dynamic 
analysis;  FEABL-5  will  be  described  in  an  ASRL  TR  in  preparation. 
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{u  . U . ...  U .}  = jth  eigenvector  (mode  shape). 
1]  2}  hj 


jth  eigenvalue  = (w_.  ) . 


jth  natural  frequency  (rad/sec) 


Total  number  of  unconstrained  degrees  of  freedom  in  the  assembled 
finite-element  model. 


Note  that  only  the  unconstrained  degrees  of  freedom  are  included  in  the 
assembled  equation  system.  If  the  constrained  degrees  of  freedom  were  also 
assembled,  the  solutions  would  include  spurious  eigenvalues  and  mode  shapes 
associated  with  the  constraints.  The  necessary  modifications  to  avoid  this 
condition  have  been  made  to  FEABL-2  by  establishing  an  internal  renumbering 
system  which  assigns  negative  numbers  to  all  constrained  degrees  of  freedom 
and  which  skips  assembly  of  the  corresponding  stiffnesses  and  masses. 

The  free-vibration  eigenvalue  problem  is  solved  by  the  Subspace  Iteration 
Method  (SIM),  generally  following  the  earlier  code  [2,3]  and  developments 
by  other  investigators  [ la] - The  SIM  has  previously  been  proposed  as  a 
computationally  efficient  method  for  analyses  in  which  only  the  lowest  few 
modes  in  the  structure  are  sought  [ 18] . Let  p be  the  actual  number  of 
eigenvalues  wanted.  Then,  generally,  the  first  P eigenvectors  are  included 
in  the  numerical  iteration,  where: 


U F f zp 

Equations  119  and  120  are  replaced  by  their  truncated  forms: 

y<„.,>  ■ IW  fei  fell 
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(123) 


However,  note  that  each  eigenvector  U_.  in  Eq.  122  still  retains  its  full  set 
of  N components. 

The  SIM  algorithm  begins  with  an  initial  guess  for  the  collection 

of  P eigenvectors.  The  eigenvectors  are  then  recomputed  in  each  major 
iteration  step  in  accordance  with  the  following  procedure.  Let  be 

the  results  for  the  eigenvectors  at  the  end  of  the  ith  iteration  step.  Then 
stiffnesses  k and  masses  m in  the  solution  subspace  are  computed  for  step 
(i+1)  by  a similarity  transformation  from  the  assembled  equations: 


Within  the  subspace,  the  auxiliary  eigenvalue  problem 


is  now  solved  by  double  Jacobi  iteration  (Subsection  7.3).  The  subspace 
matrices  k and  m are  in  general  fully  populated,  even  if  M in  the  assembled 
equation  system  is  a diagonal  (lumped-mass)  matrix.  The  collection  of 
subspace  eigenvectors  u is  a PxP  matrix,  i.e.: 


(126) 


The  matrix  operation  U u can  be  recognized  as  a vector  projection  in  which 


each  eigenvector  IP 


(i) 


is  aligned  with  its  corresponding  subspace  eic-jnvector 


u . : 
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(128 


The  remainder  of  Eq.  127  is  analogous  to  the  principal  operation  found  in 
the  simpler  Matrix  Iteration  Method  [19]  which  is  often  used  for  computation 
of  the  first  eigenvalue  only,  i.e.: 

ini ) .-l  x„  ,, 
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The  steps  outlined  in  Eqs.  124  through  127  are  repeated  until  one  of 
two  conditions  occurs.  The  iteration  will  be  stopped  after  some  maximum 
number  of  iterations  specified  by  the  user,  or  when  some  specified  number 
of  eigenvalues  n (pj<n_<P)  have  converged.  The  convergence  criterion  is 
taken  as: 
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where  the  tolerance  parameter  c.  is  specified  by  the  user.  (A  value  £=10 

will  assure  convergence  of  the  natural  frequencies  up  = to  10  V)  If, 

during  the  iteration,  some  of  the  eigenvalues  A ,A  , ...,A.  (j<n)  have 

1 z j 

already  converged,  then  the  corresponding  eigenvectors  U ,U  , ...,U.  are 

~ -L  ~ 2 — 3 

"frozen"  while  iteration  of  U_.+^,  proceeds.  The  process  is  such 

that  the  eigenvalues  computed  in  the  subspace  converge  to  the  eigenvalues  A 

of  the  assembled  equation  system.  Hence,  no  distinction  has  been  made  in 

the  notation  for  these  quantities. 

The  SIM  algorithm  outlined  by  Eqs.  124  through  127  is  not  the  most 

efficient  possible  procedure,  since  it  requires  three  coefficient  matrices 

from  the  assembled  equation  system:  K,  M,  and  K 1 in  terms  of  its  factored 

T 

triple  product  K = LDL  for  simultaneous  solution  [6].  However,  the  need 


for  the  unfactored  stiffness  matrix  K can  be  eliminated  by  the  following 
revised  algorithm.  Let  be  a collection  of  auxiliary  vectors  formed  at 

the  end  of  the  ith  iteration  step,  such  that  K = V^.  Then  the 

current  eigenvectors  are  first  calculated  from: 


(131) 


Equations  124  can  now  be  replaced  with  the  following  seris  of  computations 
which  never  require  K directly: 


(132) 


(133) 


(134) 


The  eigenvalue  problem  in  the  subspace  is  then  solved  in  the  manner  mentioned 
previously,  and  the  iteration  cycle  is  completed  with: 


(135) 


(i+1) 


It  is  now  easy  to  see  that  the  definitions  are  consistent,  i.e.  V = K U 

by  comparing  Eq.  135  with  Eq.  127. 

The  procedure  outlined  by  Eqs.  131  through  135  requires  only  the 

II  _ Ilf  rp 

assembled  mass  matrix  M and  K (in  the  form  L D L ) , and  is  thus  more 
efficient  than  the  simpler  iteration  procedure  which  was  presented  in  the 
first  part  of  this  subsection.  However,  note  that  since  both  procedures 
require  K ^ , the  SIM  is  restricted  to  analyses  of  finite-element  models  to 
which  displacement  boundary  condtions  have  been  properly  applied  to  restrain 
all  rigid-body  modes. 


7.3  Eigenvalue  Solution  in  the  Subspace 

The  eigenvalue  problem  in  the  subspace  (Eq.  125)  is  solved  by  means  of 
a two-stage  iteration  scheme  in  which  each  stage  treats  a standard  eigenvalue 


problem  with  one  coefficient  matrix.  The  eigenvalues  and  eigenvectors  of 
the  mass  matrix  m alone  are  computed  first,  and  these  quantities  are  then 
used  to  transform  Eq.  125  to  a standard  problem  [20] . Let  the  first 


eigenvalue  problem  be  represented  by: 


where 


= * y 
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= jth  eigenvalue  of  m 

After  z and  y have  been  computed,  the  mass  matrix  may  be  expressed  in  the 


Also,  powers  of  m can  be  expressed  in  simple  form,  in  particular  [20] : 
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The  validity  of  Eqs.  140  can  be  proved  in  a straightforward  manner  by 

. , . . J 1/2  1/2  -1/2  -1/2 
verifying  that  the  matrix  products  m m and  m m are  equal  to 

m and  m \ respectively.  The  proof  proceeds  by  substitution  and  carrying 

out  the  matrix  multiplications,  recognizing  that  the  collection  of  eigen- 

T T 

vectors  are  orthonormal,  i.e.  zz  =zz=  identity  matrix. 

In  order  to  transform  the  subspace  eigenvalue  problem  to  a standard 
problem,  let  v = m^^u  and  premultiply  Eq.  125  by  m to  obtain: 


-J/i  /,  -V*  if  \ -*fx  ( 'U  't.._  , 

(km  mu]  = ™ ' 2 n\  u A 


( m U k w 'L  ) V = V >* 


(142) 


Thus,  the  second  standard  problem  consists  of  forming  the  coefficient  matrix 
-1/2  -1/2 

m k m and  computing  the  eigenvalues  A and  transformed  eigenvectors  v. 

When  the  problem  has  been  solved,  the  proper  eigenvectors  u are  obtained  from 
the  definition  of  v,  i.e.: 


U = m '*  V 


(143) 


The  two  standard  eigenvalue  problems  represented  by  Eqs.  136  and  142 
are  solved  by  the  well  known  Jacobi  Iteration  Method  [19].  The  fundamental 
operation  in  the  Jacobi  method  is  a rotation  transformation  which  is  used 
to  diagonalize  the  coefficient  matrix.  Only  one  such  operation  is  required 
for  a 2x2  matrix.  Using  the  mass  matrix  m as  an  example. 
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where 


6 = - t« 


™st 


(145) 


.!  j.t  matrices,  the  procedure  is  approximate.  Successive  rotation 
■ fi.it  i vis  must  be  applied  to  reduce  the  sum  of  squares  of  the  off- 
Let  m represent  the  current  matrix  and  suppose  that  the 
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term  is  to  be  reduced  to  zero  during  the  current  step.  This  is 

accomplished  by  replacing  m with  m = R^m  R,  where: 
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The  replacement  operations  implied  by  the  product  R m R are  as  follows: 
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It  can  be  shown  [19]  that: 
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In  other  words,  the  sum  of  squares  of  the  diagonal  terms  in  m increases  by 


2(m^p)  as  a result  of  the  Jacobi  rotation  and  hence,  the  sum  of  squares  of 


the  off-diagonal  terms  must  decrease  by  that  amount.  However,  other  off- 


diagonal  terms  (m  . , m„ . ) which  may  have  been  reduced  to  zero  by  previous 
k]  x,  j 
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operations  are  made  non-zero  by  the  current  operation  which  reduces  m^. 
Thus,  the  entire  process  involves  a series  of  sweeps  through  the  coefficient 
matrix,  each  sweep  consisting  of  successive  reduction  operations  applied  to 
all  off-diagonal  terms,  until  either  a maximum  allowed  number  of  iterations 
is  reached  or  until: 


E (Current  off-diagonal  terms)  /£  (Original  off-diagonal  terms)  < £ (149) 


where  the  number  of  iterations  and  the  tolerance  parameter  £ are  specified 
by  the  user.  The  series  of  Jacobi  rotations  which  reduce  m to  a nearly 
diagonal  matrix  can  also  be  used  to  compute  the  eigenvectors  directly,  as 
shown  by  the  following  identity.  Suppose  that  J such  rotations  were  required 
in  the  order  R , R , ...,  R . Then  (see  Eq.  136): 


(150) 


The  product  for  z is  accumulated  as  the  individual  Jacobi  rotations  are 
applied. 

The  Jacobi  method  outlined  in  Eqs.  144  through  150  is  applied  in 

exactly  the  same  manner  to  the  second  standard  problem,  in  which  the 
-1/2  -1/2 

eigenvalues  of  m km  are  sought.  The  entire  procedure  has  been 
programmed  as  FEABL-5  subroutine  JACKM1,  which  employs  IBM  Scientific 
Subroutine  Package  subroutine  EIGEN  for  the  Jacobi  iteration  computations. 
Double-precision  arithmetic  is  used  within  the  subspace  to  obtain  the 
necessary  accuracy  in  the  computations  for  the  rotation  angles  0 (Eq.  146) . 
Single-precision  calculations  were  tried,  but  were  found  to  make  the  con- 
vergence process  slower,  in  agreement  with  the  conclusions  reached  by  other 
investigators  [19].  The  inefficiencies  of  double-precision  arithmetic  and 
fully  populated  matrices  are  acceptable  within  the  subspace,  which  will 
usually  be  of  the  order  20x20  or  less.  This  has  no  effect  upon  the  global 
computations,  since  K and  M are  kept  as  single-precision  band-matrices. 


7.4  Transient  Response  Analysis 

"'he  Modal  Superposition  Method  (MSM)  has  been  chosen  for  transient 
response  analysis  in  order  to  avoid  the  time-step  limitations  which  are 
imposed  by  stability  criteria,  frequency  distortion,  and  false  damping 
associated  with  finite— difference  time  operators  employed  in  direct  time- 
integration  schemes  [21,22].  In  the  direct  time- integration  approach, 
the  assembled  transient  equation  system 


M ‘j.  (t)  i K fit)  = £>(+) 

is  replaced  by  a finite-difference  equivalent,  e.g.: 

(it)*-  ^ 2li  + jri-i  ) + K fi  = Qi 


(151) 


(152) * 


where  the  subscripts  i-1,  i,  i+1  denote  quantities  evaluated  at  the  times 
+ *'=i'^+jL=i-^+At.  The  displacement  solution  q^+^  can 

then  be  estimated  from  current  and  prior  information: 


i»vi = zh~  h-' + + 

(Equation  152  applies  to  the  second  and  all  succeeding  time  steps.  For 

the  first  step,  q^  is  calculated  from  q^,  and  q^.)  Direct  time-integration 

has  the  advantage  of  not  requiring  the  assembled  equation  eigenvalue  solutions 

^or  transient  response  calculations.  However,  the  integration  stability 

boundary  does  require  knowledge  of  the  highest  frequency  w in  the  finite- 

N 

element  model,  since  At  < 2/u^  is  the  stability  criterion  [23]  for  the  central 
difference  operator.  Thus,  even  direct  time— integration  must  be  preceded 
by  at  least  the  simple  inverse  iteration  scheme: 


(153) 


U(ia)  = m'‘ku;  fi/aj 


(154) 


(Compare  Eq.  154  with  Eq.  129.) 

Since  some  knowledge  about  the  frequencies  and  mode  shapes  of  a structure 
is  usually  sought  anyway,  the  MSM  offers,  potentially,  a much  more  efficient 

* 

The  central  difference  operator  is  used  here  as  a simple  illustrative 
example . 


solution  scheme.  In  this  approach,  the  assembled  transient  system  (Eq.  151) 
is  first  decoupled  by  means  of  the  system  eigenvectors.  The  transient  dis- 
placement vector  q (t)  can  be  represented  approximately  as  a linear  combination 
of  the  first  few  eigenvectors: 

I3 

alt)  =■  I U.  e<:(t)  = U ai(-t)  (155) 

I >1  ~J  J 

where  a(t)  = {a,  a ...  a } is  the  modal  amplitude  vector.  It  is  assumed 

1 2 p 

that  only  the  p modes  actually  sought  in  the  eigenvalue  analysis  (see 

Subsection  7.2)  are  retained  in  the  modal  solution.  This  illustrates 

immediately  the  inherent  disadvantage  of  the  MSM,  that  the  quality  of  the 

solution  degrades  rapidly  for  transient  load-histories  whose  power  content 

for  frequencies  oj  > to  is  significant.  If  Eq.  155  is  now  substituted  in 

P T 

Eq.  151  and  the  result  is  premultiplied  by  U the  following  decoupled 
equation  system  is  obtained: 


UT  M (u  a)  + uT  K (U*  ' - Vi 


+ k oi  = 4 1 


(156) 


where 
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< r 


m 
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T 

= U MU 


= U K U 


f(t)  = U Q(t) 


= subspace  mass  matrix  (diagonal) 

= subspace  stiffness  matrix  (diagonal) 
= generalized  subspace  force  vector 


The  subspace  mass  and  stiffness  matrices  are  diagonal  in  the  present  case 
because  they  result  from  computations  with  a converged  set  of  eigenvectors. 

For  practical  purposes,  the  diagonal  matrices  m and  k are  obtained  simultaneously 
with  U at  the  conclusion  of  the  SIM  algorithm  (Subsection  7.2).  Note  that 
while  m and  k are  in  general  fully  populated  matrices  during  execution  of 
the  SIM,  they  become  diagonalized  as  a part  of  the  convergence  process. 

Equation  156  may  now  be  rewritten  as  the  collection  of  scalar  equations: 


\-Vf)  + kJ  *j(t)  * fiM  Qxbz>  ■■■> 


(157) 
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k./m.  (and  dropping  the  subscript  j without  loss  of  generality) 


(158) 


Damping  should  be  added  to  the  transient  response  equation  in  order  to 
permit  realistic  solutions  to  be  calculated.  This  is  done  most  conveniently 
in  the  MSM  by  directly  assuming  nondimensional  modal  damping  factors: 


(159) 


rather  than  assuming  the  physical  damping  coefficients  c 
the  typical  modal  amplitude  response  is  thus  replaced  by 


(160) 


A general  solution  of  Eq.  160  can  be  derived  by  analytical  integration 
if  a specific  assumption  is  made  about  the  generalized  force  history  f (t) . 

A piece-wise  linear  representation  for  f (t)  is  a simple  approximation  which 
is  also  useful,  since  it  is  unlikely  that  the  force  time-history  is  known 
to  an  accuracy  which  would  justify  higher-order  representation.  Therefore, 
during  the  ith  time  step  the  generalized  force  is  taken  as: 


(161) 


, = f(t.  , ) and  f.  = f(t.).  With  Eq.  161  substituted  for  f(t)  in 

1 l- 1 li 

the  following  general  solution  may  be  derived: 


(162) 


(163) 


where 
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= a(t. 
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amplitude  at  beginning  of  time  step 
amplitude  velocity  at  beginning  of  time  step 
damped  natural  frequency 


and  where 


c = 


I.f.  ' 
li  J hl-  > 


(164) 


Equations  162  and  163  give  the  modal  amplitude  on  and  its  velocity  on 

at  the  end  of  the  time  step  t=t^,  in  terms  of  the  associated  generalized 

force  history  (f.  . f.)  prior  information  (a.  ,,  a . ) and  system  properties 

l-l  i l-l  i-I 

(£,0),k).  It  should  be  noted  that  this  solution  applies  only  to  underdamped 
systems,  i.e.  L,  < 1. 

The  MSM  solution  scheme  is  quite  straightforward.  Vectors  of  initial 
modal  amplitudes  and  modal  amplitude  velocities  rx^,  CXQ  are  prescribed  and 
the  set  of  discrete  force  vectors: 

[{0  fL  l2  J = yTf<2^1  •••}  (165) 


is  computed.  Equations  162  and  163  are  then  employed  for  each  time  step  in 
succession,  and  the  physical  displacement  solutions 


f 


(166) 


are  calculated  during  the  transient  solution.  The  scheme  presented  in 
Eqs . 162  and  163  has  the  advantage  that  the  size  of  each  time  step  may  be 
arbitrarily  chosen.  The  most  economical  choice  of  times  t^ , t^,  ...  is 
such  that  each  interval  is  as  long  as  possible  while  still  maintaining  a 
reasonably  accurate  piece-wise  linear  approximation  for  the  force  history. 
Specific  choices  are,  of  course,  dependent  upon  the  specific  problem  under 
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analysis.  Tne  transient  solution  algorithm  lias  been  programmed  as  FEABL-5 
subroutines  TSTEP  and  ELAPSE. 

7.5  Software  Verification  and  Performance  Tests 


Several  test  analyses  were  run,  in  which  the  computed  results  were 
compared  with  independent  analytical  solutions,  both  to  verify  the  basic 
software  and  to  assess  the  performance  of  the  SIM.  The  verification  tests 
included  validations  of  subroutine  MLP3M,  the  element  mass  matrix  generator 
which  is  the  companion  to  element  stiffness  matrix  subroutine  MLP3K.  Results 
of  the  tests  are  presented  here  approximately  in  chronological  order,  i.e. 
following  the  actual  development  work. 

7.5.1  Verification  of  Jacobi  Iteration  Method 

The  accuracy  of  the  Jacobi  Iteration  Method  (subroutines  JACKM1 
and  EIGEN,  Subsection  7.3)  was  assessed  in  the  first  series  of  tests  by 
computing  the  first  few  natural  frequencies  for  both  cantilever  and  unre- 
strained slender  beams.  The  purpose  of  this  series  of  tests  was  two-fold: 
first  to  reassess  the  sensitivity  of  the  Jacobi  method  to  arithmetic  pre- 
cision [19]  and  second,  to  verify  the  code  in  subroutine  EIGEN.  Some  doubt 
had  been  cast  upon  the  validity  of  parts  of  the  IBM  Scientific  Subroutine 
Package  because  it  is  no  longer  supported  [24].  However,  the  verification 
tests  did  demonstrate  the  validity  of  subroutine  EIGEN,  as  will  be  seen 
below. 

Independent  analytical  solutions  of  the  slender-beam  vibration 
problem  [20]  give  the  following  results  for  the  first  few  natural  frequencies 
of  a cantilever  beam  with  uniform  section  properties: 

<Vi  * (1.17  5)*  -ft 


= (i.mf  vll'7wL4 
= (7.8'r*r)i'  /e  i /m  \F 


where  the  frequencies  are  in  units  of  rad/sec,  and  where: 
E -=  material  Young's  modulus  (psi) 

4 

I = cross  section  bending  inertia  (in  ) 

2 2 

m = mass  per  unit  length  (lb  sec  /in  ) 

L = total  length  of  the  beam  (in.) 


Similarly,  the  solutions  for  an  unrestrained  beam  are: 


- — j 


= U)2  - 0 

=■  (l.S'Obtr)2'  •/ E 3 /rv)  L*  (168) 
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The  first  two  solutions  in  Eq.  168  correspond  to  the  rigid-body  modes  of 
vertical  translation  and  rotation  of  the  beam  about  its  center. 

A cubic-interpolation  assumed-displacement  beam  element  was 
programmed  to  provide  element  software  to  implement  the  tests.  The  beam 
element  (Fig.  15)  is  of  length  SL,  such  that  the  total  beam  length  L is  a 
multiple  of  the  element  length.  The  element  stiffness  matrix,  which  is 
easily  derived  by  substituting  the  displacement  interpolation  into  the 
strain-energy  equivalence: 
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is  given  by: 
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The  representation  is  "exact",  in  the  sense  that  the  engineering  beam  theory 
solutions  for  static  deflections  of  a cantilever  under  a tip  load  or  tip 
moment  (for  example)  can  be  obtained  exactly  with  models  in  which  the  entire 
beam  is  represented  by  a single  element.  In  other  words,  the  beam  element 
introduces  no  convergence  error  due  to  discretization  of  a continuum.  However, 
some  convergence  error  in  the  higher  eigenvalues  can  still  be  expected,  even 
with  several  elements  modelling  the  beam,  as  will  be  explained  subsequently. 

The  element  consistent  mass  matrix  is  derived  by  substituting  the  interpolation 
function  into  the  kinetic  energy  equivalence: 
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from  which  there  results 


(172) 


Analytical  and  computed  results  are  compared  for  both  the 
cantilever  and  unrestrained  beams  in  Table  4.  Computed  values  are  shown 
for  both  single-  and  double-precision  arithmetic.  It  is  evident  that  the 
arithmetic  precision  has  not  affected  the  results  except*  in  the  case  of 
the  rigid-body  modes  of  the  unrestrained  beam.  Therefore,  single-precision 
arithmetic  could  be  used  when  subroutines  JACKM1  and  EIGEN  are  coupled  with 
the  SIM  software,  since  the  SIM  can  be  applied  only  to  structures  which 
are  restrained  against  any  rigid-body  motion.  However,  the  results  for  the 
unrestrained  beam  do  indicate  that  accuracy  problems  might  occur  for  struc- 
tures which  possess  a very  low  first  natural  frequency,  or  for  structures 
which  possess  one  or  more  closely  grouped  frequencies.  The  arithmetic  in 
subroutines  JACKM1  and  EIGEN  has  been  left  in  double-precision  for  the 
above  reason  and  to  permit  their  use  as  an  independent  module  for  eigenvalue 
analyses  of  unrestrained  structures.  The  apparent  "improvement"  in  the 
single-precision  results,  in  comparison  with  the  earlier  investigation  which 
gave  rise  to  the  warning  about  precision  [19],  is  most  probably  due  to  the 
transition  from  12-  and  16-BIT  hardware  in  the  late  1950' s to  today's 
32-BIT  hardware. 

The  errors  in  both  computed  solutions  for  the  higher  modes 
serve  to  illustrate  another  form  of  convergence  error  which  is  specifically 
associated  with  eigenvalue  analysis.  As  higher  natural  frequencies  are 
sought  for  the  beams,  the  corresponding  exact  mode  shapes  exhibit  higher 


It  was  also  noted  that  more  iterations  were  required  in  single-precision 
to  achieve  a given  error  tolerance  e. 


spatial  frequency  distributions  of  transverse  displacement  as  a function  of 
distance  along  the  beam.  The  ability  of  the  finite-element  model  to  reproduce 
mode  shapes  and  frequencies  is  thus  limited,  through  the  element  displacement 
interpolation  function,  by  the  number  of  elements  in  the  model.  This  point 
is  further  illustrated  by  the  results  in  Table  5,  which  compares  two  solutions 
for  the  cantilever  beam.  Here,  the  arithmetic  precision  has  been  kept  constant 
but  one  model  contains  only  12  degrees  of  freedom  (5  elements)  while  the  other 
contains  50  degrees  of  freedom  (24  elements) . The  improvements  are  apparent 
in  the  third  and  higher  frequencies. 


The  SIM  software  was  first  verified  by  conducting  another  series 
of  tests  in  which  SIM  solutions  of  the  cantilever  beam  problem  were  compared 
with  analytical  results  outlined  in  Subsection  7.5.1.  A 50  degree-of- freedom 
model  was  used  to  generate  the  assembled  equation  system.  The  size  of  the 
subspace  was  set  at  P=12,  while  convergence  was  sought  for  p=6  eigenvalues. 

The  results  presented  in  Subsection  7.5.1  were  reproduced. 

In  a second  series  of  tests,  the  SIM  software  was  coupled  with 
element  MLP3K  and  its  mass  matrix  MLP3M  (Section  4)  to  provide  further  verifi- 
cation while  simultaneously  validating  subroutine  MLP3M.  The  single-layer, 
isotropic,  simply  supported  square  plate  illustrated  in  Fig.  16  was  chosen 
as  a test  problem.  One  quadrant  of  the  plate,  shown  in  the  lower  part  of 
the  figure,  was  isolated  for  analysis  by  applying  symmetry  boundary  conditions 
along  the  centerlines.  Eigenvalue  analysis  of  a model  of  this  type  will 
compute  only  those  frequencies  corresponding  to  mode  shapes  which  are 
symmetric  in  both  the  X and  Y directions.  Independent  analytical  solutions 
for  the  symmetric  frequencies  were  obtained  from  a Fourier  analysis  [25] 
for  comparison.  A 3x3  mesh  is  shown  in  Fig.  16,  but  5x5,  7x7  and  9x9  meshes 
were  also  analyzed  to  obtain  information  about  solution  convergence  rates. 

In  the  first  test,  the  subspace  size  was  set  at  P=9,  the  error 

-2 

tolerance  parameter  was  set  at  E = 10  , and  convergence  was  sought  for 

p=6  eigenvalues.  The  results  are  summarized  in  terms  of  percent  errors: 


(173) 


Computed  Exact 


and  are  shown  in  Table  6 and  Fig.  17.  Analyses  were  run  using  both  hybrid- 
rational  and  lumped  element  mass  matrices.  Surprisingly,  the  lumped-mass 
matrix  appears  to  converge  as  fast  as  or  faster  than  the  hybrid-rational 
mass  matrix  for  the  higher  frequencies,  while  the  hybrid-rational  mass 
matrix  converges  faster  for  the  first  natural  frequency.  This  result  is 
opposite  to  what  one  would  expect  from  a comparison  of  a lumped-mass  matrix 
with  a consistent-mass  matrix  for  an  assumed-displacement  element.  In 
either  case,  the  convergence  is  quite  rapid  after  enough  detail  has  been 
placed  in  the  finite-element  model. 

A more  interesting  conclusion  may  be  drawn  from  the  erroneous 

predictions  associated  with  the  combinations  of  higher  frequency  and 

coarser  mesh.  Here  reappears  the  problem  of  mesh  convergence  error  which 

was  mentioned  in  Subsection  7.5.1.  The  results  are  dramatic  in  the  present 

case  because  of  the  linear  interpolation  which  is  used  in  element  MLP3K  for 

the  plate's  transverse  displacement,  w.  The  convergence  effect  is  illustrated 

schematically  in  Fig.  18,  which  compares  the  spatial  distributions  of  the 

first  few  natural  modes  with  the  abilities  of  the  finite-element  models  to 

mimic  these  shapes.  Fcr  example,  the  3x3  mesh  is  able  to  reproduce  the 

first  and  third  harmonics  without  spatial  distortion,  but  can  not  do  so  for 

the  fifth  harmonic.  Thus  (considering  the  lumped-mass  results  in  Table  6), 

low  errors  appear  for  f , f^  and  f ^ . The  next  eigenvalue,  f^,  should 

also  have  a low  error,  but  it  is  apparently  "infected"  by  its  next  neighbors 

f._  and  f . A similar  argument  can  be  made  about  the  5x5  mesh,  which  is 
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capable  of  reproducing  harmonics  up  to  the  fifth,  but  which  severely  distorts 
the  seventh  harmonic. 

A second  test  was  conducted  to  investigate  the  effect  of  subspace 

size,  P,  on  the  rates  of  convergence  of  the  six  lowest  eigenvalues.  The 

finite-element  mesh  was  fixed  at  9x9  for  this  test,  in  order  to  eliminate 

as  nearly  as  possible  the  effects  of  modelling  convergence  error.  As  can 

be  seen  in  Table  6,  the  worst  errors  with  a 9x9  mesh  and  P=9  were  7.08 

percent  with  the  hybrid-rational  mass  matrix  and  2.29  percent  with  the 

lumped-mass  matrix.  Subspace  sizes  P=6,8,10  and  12  were  investigated  in 

the  second  test,  for  comparison  with  the  case  P=9.  The  tolerance  parameter 

-2 

was  kept  fixed  at  C = 10 


According  to  the  mathematical  theory  of  eigenvalue  analysis  [18] , 
the  principal  reason  for  the  potential  efficiency  of  the  SIM  as  a numerical 
method  is  that  the  SIM  algorithm  causes  eigenvalues  A to  converge  at  rates 
proportional  to  A /Ap^,  while  the  rates  for  other  numerical  methods  may 


only  be  proportional  to  A 
sought  and  if: 


Thus,  for  example,  if  6 eigenvalues  are 


then  the  choice  P=12  should  be  superior  to  P=10  because  A /A ^ ^ < A^A^ 
(j=l,2,...,  6).  This  situation  occurs  for  the  present  plate  analysis,  as 
illustrated  by  the  analytical  eigenvalue  solutions  in  the  second  column  of 
Table  7.  The  remainder  of  the  table  gives  the  nondimens ional  convergence 
rate  factors  A^/A^^  (j=l»2,...,  6)  for  the  various  choices  of  P which  were 
investigated.  Theoretically,  the  error  in  a computed  eigenvalue  at  any 
iteration  step  may  be  expressed  as: 


X 

Uj/Ap+i)  - (Rj) 


where  I represents  the  cumulative  number  of  iterations  which  have  been 
completed.  Taking  the  logarithm  of  Eq.  175  and  recognizing  that 

I - Elapsed  CPU  Time 

then  leads  to  the  relation: 


log 10  (Error) 

— — r — = Constant  x log.,-(R.) 

(CPU  Time)  10  g 


The  results  of  the  test  are  summarized  in  Table  8 in  terms  of 
percent  errors  in  the  predictions  for  the  first  6 eigenvalues.  The  CPU  times 
varied  approximately  from  0.38  to  0.66  minute  per  case,  as  indicated  in  the 
bottom  row  of  the  table.  Shown  below  each  error  is  a figure  in  parentheses 
equal  to  the  quantity  log  (Error)/ (CPU  Time)  in  arbitrary  units.  A cross-plot 
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of  the  data  in  Tables  7 and  8 should  lead  to  a linear  relationship  for  each 
type  of  element  mass  matrix  on  semi-logarithmic  paper  if  Eq.  176  is  valid. 
Such  a plot  is  given  in  Pig.  19,  and  it  indicates  that  the  simple  convergence 
theory  outlined  above  is  not  obeyed  by  the  present  case.  The  probable  cause 
of  the  discrepancy  is  that  the  theoretical  convergence  rates  are  asymptotic 
values  [18],  i.e.  rates  which  are  realized  near  the  end  of  the  entire  SIM 
iteration  process,  while  a significant  portion  of  the  CPU  time  is  expended 
upon  iterations  for  which  the  convergence  rates  are  not  asymptotic.  However, 
the  actual  CPU  times  of  less  than  one  minute  for  eigenvalue  analyses  of 
500  degree-of- freedom  models  still  demonstrate  that  the  SIM  is  an  efficient 
computing  algorithm. 

7.5.3  Verification  of  Transient  Response  Analysis 


The  transient  response  software  (subroutines  TSTEP  and  ELAPSE, 
and  Eqs . 162  and  163;  see  Subsection  7.4)  was  verified  by  four  test  analyses 
of  the  3x3-mesh  finite-element  model  of  one  quadrant,  of  a simply  supported 
square  plate  having  the  following  properties: 


E = 10  psi 
Single  layer. 


V = 0.3 


h = 0.01  inch 


Mass  density, 


0.1  lb  sec  /inch 


Edge  dimension, 


L = 10  inches 


The  first  two  tests  attempted  to  reproduce  the  correct  response 
to  a steady-state  sinusoidal  forcing  function  given  by  the  pressure  distribu 
tion: 

\ , I it  X \ /ir  Y \ . / \ 

y,t)  = fc  l ~j~/  Sir)  (.  j_  ) J.n  (Jit/  (178) 

with  p =0.5  psi.  The  above  distribution  is  designed  to  excite  the  first 
o 

natural  mode  of  the  plate,  for  which  the  transverse  displacement  shape 
function  is  given  by 


/'  TA  \ it  v ' 

Kl  S t n ( £7  n n ( -j 
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The  corresponding  values  for  the  generalized  force  and  the  generalized  mass 
are  then: 

rLr  L 

Jo  J ?IX,V,0  ^(X,y)  JXclY  --  — ^aL  si^(SLt)  (18C 


p L p L st  i ^ 

- j J oh  (xi  Y)  JX^Y  ~ ^ 


The  steady-state  solution  of  Eq.  160  gives  the  following  result  for  the 
modal  amplitude  a^,  which  is  also  numerically  equal  to  the  amplitude  of 
the  center  deflection  of  the  plate: 

y = fl  /^l  

[u-  siY'/ol  )L  + (zja/^i)2  J1/L 

where,  for  the  properties  given  by  Eqs.  177,  the  first  natural  frequency  is: 

- 6-11  rad/sec  (182 

The  first  test  was  run  at  = 0.5u^  and  £ = 0.1,  for  which 
= 10.88  in.  according  to  Eq.  182.  In  this  case,  the  computed  center 
deflection  amplitude  was  found  to  be  10.33  in.,  giving  an  error  of  5.05 
percent.  The  second  test  was  run  at  Q = and  £ = 0.05,  for  which 
= 82.4  in.  according  to  Eq.  182.  In  this  case,  the  computed  center 
deflection  amplitude  was  found  to  be  80.9  in.,  giving  an  error  of  1.82 
percent.  Plots  of  both  runs  indicated  that  the  center  deflection  amplitude 
built  up  to  its  steady-state  value,  as  should  be  expected.  Figure  20 
illustrates  the  computed  center  deflection  time  history  from  the  second 


The  third  and  fourth  tests  were  run  to  verify  the  initial- 
condition  portions  of  the  software.  With  the  damping  factor  kept  constant 
at  £ = 0.1,  the  plate  was  subjected  to  the  initial  conditions: 

w(X,Y,  o)  = o.s  wt  (xy)  w(X,Yl0)*O  (IE 


J 
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in  the  third  test,  and 


w 


=o 


w (X,Y,  o)  = 0.5T  STj  (*, 


(185) 


in  the  fourth  test,  with  the  forcing  function  p(X,Y,t)  = 0 in  both  tests, 
in  order  to  simulate  damped  free  vibration.  Time  histories  of  the  computed 
center  deflections  from  these  runs  are  illustrated  in  Figs.  21  and  22, 
respectively.  The  behavior  has  the  proper  appearance,  and  the  calculations 
summarized  in  the  figures  show  that  the  input  damping  factor  £ = 0.1  has 
been  reproduced  by  the  computed  results  with  about  3 to  4 percent  error. 

The  solutions  are  probably  somewhat  better  than  these  results  indicate, 
since  the  estimations  of  £ have  been  based  upon  roughly  sketched  decay 
envelopes. 
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SECTION  8 


I ' 


APPLICATIONS  PROGRAMS 


8.1  Introduction 

Several  applications  programs  were  prepared  during  the  investigation 
in  order  to  verify  software  integration  and  to  address  some  specific  problems 
of  current  interest.  Each  applications  program  is  written  in  the  form  of 
a primary  subroutine,  which  controls  mesh  generation  and  execution  of  the 
analysis,  and  a dummy  main  program  which  simply  calls  the  primary  subroutine. 
The  purpose  of  the  dummy  main  program  is  to  permit  the  re-dimensioning  of 
certain  FORTRAN  vectors  and  arrays,  when  required,  without  the  need  for 
recompilation  of  the  more  extensive  code  in  the  primary  subroutine. 

Each  applications  program  is  tailored  to  employ  only  those  general 
software  modules  specifically  required  for  the  particular  analysis  task. 

Block  diagrams  of  the  individual  general  modules  are  illustrated  in  Fig.  23. 

Those  application  programs  which  are  expected  to  be  useful  in  the 
future  are  described  in  this  section.  The  appendices  contain  listings  of 
the  dummy  main  and  primary  subroutine  for  each  program  described  below. 
Modular  block  diagrams  in  each  subsection  summarize  the  general  software 
required  by  each  applications  program. 

8.2  Subroutine  AMMRC  and  AMMRC2 

One  of  the  objectives  of  the  present  investigation  was  to  assess  the 
need  for  special  traction- free-edge  versions  of  the  hybrid  plate  elements 
for  the  purpose  of  obtaining  accurate  stress  solutions  at  free  edges  of 
structural  details.  Hybrid  elements  without  this  special  feature  are 
usually  capable  of  reproducing  the  traction-free  condition  in  geometrically 
simple  situations,  e.g.,  zero  bending  and  shear  stresses  at  the  edges  of  a 
simply  supported  rectangular  plate  subjected  to  transverse  pressure  load- 
ing [2,3].  The  situation  is  less  clear  for  more  complex  situations,  e.g. 
the  bending  of  a rectangular  plate  which  contains  a circular  cutout. 
Therefore,  two  applications  programs  were  prepared  to  provide  analyses  of 
this  situation.  Subroutine  AMMRC  was  intended  to  model  the  structure  with 
thick-plate  elements,  using  either  regular  or  traction-free  elements  along 
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the  free  edges.  This  program  was  discontinued  because  of  the  problems 
encountered  in  computing  efficiency  and  derivation  of  the  special  traction- 
free  stress  distribution  (see  Section  3) . Subroutine  AMMRC2  (listed  in 
Appendix  A)  models  the  structure  with  moderately  thick  quadrilateral  plate 
elements  (MLP3K) , for  which  the  traction-free  formulation  of  assumed  stresses 
is  not  possible  (see  Section  4) . 

A brief  literature  search  failed  to  provide  any  independent  analytical 
solutions  for  this  problem.  Therefore,  the  finite-element  model  was  designed 
to  simulate  structure  and  loading  which  could  easily  be  reproduced  in  the 
laboratory,  in  order  to  obtain  experimental  data  for  comparison  with  the 
computed  stresses.  The  configuration,  illustrated  in  Fig.  24,  is  a 
rectangular  plate  with  a symmetrically  placed  circular  hole.  The  plate  is 
subjected  to  four-point  bending  by  four  knife  edges,  so  that  the  bending 
moment  is  constant  in  the  test  section.  A quadrant  of  the  plate  is  isolated 
for  analysis,  as  shown  in  the  lower  part  of  the  figure.  This  model  is 
restricted  to  analyses  of  isotropic  plates  and  0°/90°  balanced  crossplies 
because  of  its  symmetry  boundary  conditions,  which  are  not  valid  for  other 
laminates  possessing  bending- twisting  coupling.  The  dimensions  A,  A^,  A^, 

B,  and  R constitute  the  geometrical  input,  from  which  the  finite-element 
mesh  is  automatically  generated.  Scale  mesh  plans  from  typical  runs  are 
illustrated  in  Fig.  25. 

The  mesh-generation  algorithm  is  summarized  schematically  in  Fig.  26, 
which  shows  the  model  divided  into  four  regions.  The  first  region  is  sized 
to  have  a square  outer  boundary,  so  that  a 45°  ray  subdivides  the  region 
into  two  equal  areas.  There  are  always  8 elements  placed  circumferentially 
around  the  quarter-circle  boundary,  a condition  which  determines  that  the 
other  regions  always  have  four  rows  of  elements  from  the  plate  centerline 
to  the  lateral  free  edge.  The  edge  length  of  a typical  element  along  the 
quarter-circle  boundary  is  rr/16,  and  the  number  of  elements  allowed  radially 
(NER)  in  the  first  region  is  computed  to  make  the  aspect  ratios  of  the 
boundary  elements  as  close  to  unity  as  possible*  (see  Fig.  25) . The  second. 


It  is  sometimes  necessary  to  reprogram  this  part  of  the  code  to  increase 
the  aspect  ratios  of  these  elements,  in  order  to  reduce  aspect  ratios  in 
other  parts  of  the  mesh.  For  example,  an  aspect  ratio  of  2 was  used  for 
the  comparison  with  experimental  results  discussed  at  the  end  of  this 
subsection. 


I 


; 


third  and  fourth  regions  are  filled  with  rectangular  meshes,  with  the 
number  of  elements  between  the  bounding  lines  of  each  region  computed  to 
make  the  element  aspect  ratios  as  close  as  possible  to  being  between  1 and 
2.  A few  of  the  key  element  and  node  numbers  are  given  in  Fig.  26.  Also, 
the  element  numbering  sequences  are  represented  by  arrows,  and  the  DO  loop 
in  the  code  which  generates  each  region  is  noted  to  provide  a reference  to 
the  program  listing. 

The  following  user  actions  are  required  to  execute  the  analysis.  First, 
vector  and  array  dimensions  must  be  properly  defined  in  the  dummy  main  pro- 
gram. A dimension  of  50,000  words  is  suggested  for  the  FEABL-2  data  vector 
names  (RE, IN)  as  a starting  point.  Whatever  dimension  is  chosen  must  also 
be  defined  in  the  DATA  statement: 

DATA  LENGTH,  NLY, NLY1/  xxx,  y,  z/ 


where 


xxx  = length  of  the  data  vector 

y = total  number  of  distinct  layers  in  the  laminated 
plate . 
z = y+1 


Also,  the  auxiliary  arrays  and  vectors  must  be  dimensioned  as  follows  to 
agree  with  the  DATA  statement: 

DIMENSION  XEL2 (y,6) , XR(y),  CMC(y,3,3),  HI(z),  H (y) 


Second,  job  control  instructions  external  to  the  FORTRAN  code  must  be 
provided  to  establish  two  temporary  sequential-access  datasets  on  system 
disk  or  tape  storage.  The  first  dataset  must  be  allocated  10  single- 
precision words  per  record  and  must  be  assigned  FORTRAN  unit  number*  20. 
The  second  data  set  must  be  allocated  352+9y  words  per  record  and  must  be 
assigned  FORTRAN  unit  number  21.  Each  data  set  must  be  allocated  a number 


The  FORTRAN  unit  number  (sometimes  referred  to  as  a hardware  device  code) 
is  used  to  control  input/output  operations.  The  unit  numbers  5 and  6 are 
commonly  assigned  to  the  card  reader  and  line  printer,  respectively,  on 
IBM  and  CDC  systems. 
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of  records  sufficient  to  provide  one  record  per  plate  element.  Job  control 
instructions  must  also  be  provided  to  access  the  pre-compiled  codes,  or  the 
source  decks  must  be  provided  for  the  general  software  modules  shown  in 
Fig.  27. 

Third,  the  user  must  prepare  input  data  cards  in  accordance  with  the 
conventions  given  in  Table  9.  The  geometrical  data  must  satisfy  the  restric- 
tions: 


A-A,  > B > R A . > A 1 (186) 

The  program  automatically  places  the  neutral  axis  at  the  geometrical  midplane 
of  the  plate.  If  a multilayer  laminate  is  to  be  analyzed,  the  laminate  must 
be  of  balanced  construction  (distribution  of  layer  angles  and  thicknesses 
symmetric  about  the  mid-plane)  because  of  the  symmetry  boundary  conditions 
imposed  in  the  analysis.  The  program  places  a line  load  of  1 lb/inch  along 
the  knife  edges. 

Two  preliminary  analyses  of  single-layer  isotropic  aluminum  plates  were 
conducted  with  subroutine  AMMRC2  to  assess  the  ability  of  element  MLP3K  to 
reproduce  stresses  and  strains  properly  near  the  cutout.  Stress  results  are 
summarized  in  Fig.  28  for  the  case  of  a 1-inch-thick  x 16-inch  x 21-inch 
plate  with  an  8- inch-diameter  hole  and  with  the  knife  edges  placed  3 inches 
apart.  Thus,  the  applied  line  bending  moment  per  unit  plate  width  is: 

M = 3.0  in-lb/in. 

in  the  gage  section,  and  the  maximum  bending  stress  to  be  expected  at 
locations  remote  from  the  hole  is  given  by 

<fK-  C M / H Z - 18  ps(  (187) 

The  stress  contours  shown  in  Fig.  28  demonstrate  that  Eq.  187  is  obeyed 

remote  from  the  hole,  and  also  that  the  bending  stress  decreases  linearly 

to  zero  between  the  knife  edges,  as  it  should.  More  significantly,  the 

bending  stress  approaches  zero  in  the  region  near  the  hole  and  the  X-axis, 

i.e.,  reproducing  the  traction-free  condition  at  the  hole.  Along  the  Y-axis, 

a stress  concentration  appears  in  0^,  as  should  be  expected.  An  extrapolation 

of  these  results  to  the  edge  of  the  hole  is  compared  with  independent  experimental 


results  obtained  from  a handbook  [26]  in  Pig.  29.  Although  the  handbook  data 
do  not  cover  the  range  of  parameters  used  in  the  present  analysis,  it  is 
apparent  that  the  finite-element  results  are  somewhat  stiff,  i.e.  the 
computed  stresses  are  lower  than  the  probable  experimental  values. 

An  experiment  was  run  by  AMMRC  personnel  in  order  to  obtain  better 
confirmation.  The  test  specimen  was  an  aluminum  plate  with  the  following 
properties  and  dimensions: 

E = 107  psi  and  V = 0.3 
Length  (2A)  = 10  inches 

Width  (2B)  = 4 inches 

Hole  Diameter  (2R)  = 1 inch 
Thickness  = 0.26  inch 

The  computed  and  measured  strains  along  the  Y axis  are  shown  in  Fig.  30. 

8.3  Subroutines  AMMRC3  and  AMMRC4 

One  of  the  verification  tasks  which  was  defined  during  the  investigation 
involves  analysis  of  laminated  conical  shells  subjected  to  axial  compression, 
shear  and  bending,  as  shown  schematically  in  Fig.  31.  This  analysis  provides 
a test  of  the  ability  of  the  flat-plate  MLP3K(Q)  elements  to  model  curved 
shells.  Subroutines  AMMRC3  (Appendix  B)  and  AMMRC 4 (Appendix  C)  were  pro- 
grammed for  this  purpose. 

The  lower  part  of  Fig.  31  summarizes  the  conventions  which  were  adopted 
for  the  automatic  mesh-generation  scheme  in  subroutine  AMMRC 3 . A half-model 
of  the  shell  is  used,  since  the  applied  loads  and  shell  geometry  are  symmetric 
with  respect  to  the  XZ  plane.  The  half-model  is  a faceted  surface  consisting 
of  a number  of  equi-angular  gores;  the  lateral  edges  of  the  gores  are  coin- 
cident with  generators  in  the  true  cone  surface.  Each  gore  is  then  subdivided 
into  a number  of  MLP3K  elements. 

Several  systems  of  axes  must  be  dealt  with  in  order  to  generate  the 
problem  data  and  assemble  the  final  equation  system.  Coordinates  of  points 
on  the  true  cone  surface  are  described  in  the  reference  axis  system  XYZ. 

For  example,  points  on  the  generators  AB,  CD  (Fig.  31)  are  described  by: 

R = Rt  (Z/l)  + Rjl-  2/l)  (188) 
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(189) 


X=  R cos$1  Y = R sin  &1  (^e  AB) 

X ~ R r «>  s Oj  Y=  R S i nt?2  ( Eo| j e C D) 

However,  the  MLP3K  elements  must  be  given  nodal  coordinates  in  the  local 

gore  plane  system  xyz  (common  for  all  gores) , and  the  computed  element 

stiffnesses  must  then  be  transformed  to  the  collection  of  systems  x^y^z^, 

X2y2Z2  before  assembly  of  the  final  equations.  The  latter  systems  are 

oriented  such  that  the  axes  z^,z^  are  locally  normal  to  the  true  cone 

surface.  Thus,  the  local  plate-element  degrees  of  freedom  u,v,w,G  ,0  are 

x y 

transformed  to  u.,v  .w,,0  ,0  along  AB  and  u„,v„,w  .0  ,6  along  CD  to 

1 1 1 X1  Y1  2 2 2 X2  y2 
permit  consistent  assembly  to  adjacent  gores. 

The  required  transformations  are  three-dimensional  rotations  governed 

by  the  geometrical  relationships  between  the  various  axis  systems.  It  is 

easy  to  show  that: 


fxyl}-  D JX  VZ] 

where  i=l,2  and  where 

^ (Rz-  ) CoS  d coS  3 /A 
~i>  ~ ! ( Cos  &L)  /Z  S‘ n 3 


h;  y.  *;]  = 5,-fX  Y Zi 


(^2.  ^1^  Sincdcos^l  / L /A 


(190) 


sin  - s i n 01  ) /l  si  n /3 
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L (sin#2- slo6i  l/zAs-'n  -}  - L (cosfi2  - cosB1 ) /zAsinfl  (#z- i?d  ) coSfl/A 
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> (191) 
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(192) 


Equations  190  are  combined  to  provide  a transformation  from  the  local  axis 
system  xyz,  in  which  the  element  stiffnesses  are  computed,  to  the  axes  x^y^z. 
for  assembly: 
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fx_y  zi  “■  U<  y i «()  ~ T-  f yL  */} 

The  corresponding  degree-of- freedom  transformation  is: 
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and  where  i-1  or  2 according  to  whether  the  node  considered  lies  on  the  edge 
AB  or  CD,  respectively  (see  Fig.  31).  Recognizing  that  element  nodes  1,  2 
are  on  edge  AB  and  nodes  3,  4 are  on  edge  CD  according  to  the  numbering 
convention  given  in  Fig.  31  and  formally  substituting  Eq.  194  in  a strain- 
energy  expression,  one  obtains  for  the  global  stiffness  matrix  of  a 
typical  element: 


I Si 

! R 
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where  k is  the  stiffness  matrix  computed  by  subroutine  MLP3K,  and  where 
the  transformation  matrices  in  Eq.  196  are  20x20  super-diagonal  matrices. 
Note  that  the  foregoing  procedure  is  programmable  directly  from  the  geo- 
metrical relationships  between  the  various  coordinate  systems  and  that  it 
selects  directly  the  best  set  of  coordinate  systems  for  assembly  of  the 
final  equations  (z^  coordinate  normal  to  the  true  cone  surface) . Also,  the 
element  matrices  k need  be  computed  only  once,  since  all  gores  are  geo- 
metrically identical  when  viewed  in  their  own  local  planes.  However,  the 

matrices  k must  be  recomputed  for  each  gore,  since  the  transformation 
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matrices  R^,  are  functions  of  the  angular  position  of  the  gore  on  the 
cone  surface.  Finally,  note  that  the  mesh  subdivision  is  done  most  con- 
veniently in  the  reference  coordinate  system  XYZ  (Eqs.  188  and  189),  but 
that  local  coordinates  are  obtained  easily  from  the  XYZ  nodal  coordinates 
by  means  of  the  second  of  Eqs.  190. 

The  applied  loads  are  treated  in  a similar  manner,  except  that  it  is 
most  convenient  to  express  the  loads  initially  as  stress  distributions  in 
the  reference  system  XYZ.  The  applied  bending  moment  is  assumed  as  an 
axial  stress  distribution  which  varies  linearly  with  X: 


(198) 


is  the  tip  cross  section  moment  of  inertia,  and  where  H is  the  total 
thickness  of  the  shell.  This  distribution  is  only  an  approximation  when 
the  shell  is  a multilayer  laminate,  but  St.  Venant's  Principle  insures  that 
the  computed  stresses  in  elements  away  from  the  tip  of  the  shell  will 
quickly  conform  to  the  proper  distribution  for  the  laminate.  Equations  197 
and  198  may  be  combined  to  give  the  following  expression  for  the  average 
axial  stress  due  to  bending  which  acts  upon  a typical  gore  bounded  by  edges 
at  angles  0 ,8  : 


Equation  199  is  statically  equivalent  to  a pair  of  nodal  forces 
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Fj  = - M(<-o5  61  + cosdz)  p/z-n  ^ • Q=  (Bz-9i  )/z. 

which  are  applied  at  the  two  tip  nodes  of  the  gore  (points  A and  D in 
Fig.  31) . In  a similar  manner,  the  nodal  forces: 

h,  = _ t-  3 / 2-  7T  Fj  = V ( SirFBj  + Sir*  flz  ) /j  /in 
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are  derived,  corresponding  respectively  to  a constant  distribution  of 

compression  stress  a and  a parabolic  distribution  of  shear  stress  o 

z ZX 

with  maximum  at  X=0.  The  typical  nodal  force  vector  is  then  given  by: 

% = K = l Fi  0 >VF2  i 


(202) 


However,  vector  Q is  aligned  with  the  reference  axes  XYZ,  and 

Si  = Qi  S S 


(203) 


must  be  computed  for  assembly  at  points  A and  D,  respectively,  where  D , 
■and  D^  are  given  by  Bq.  192.  The  details  of  these  computations , as  well 
as  the  stiffness  transformations,  can  be  followed  in  the  program  code 
(Appendix  B) . 

The  following  user  actions  are  required  to  execute  the  analysis. 
First,  vector  and  array  dimensions  must  be  properly  defined  in  the  dummy 
main  program.  The  FEABL-2  data  vector  (variables  RE,  IN)  should  be  given 
the  dimension: 


xxx  = (zSblg  f 45 -'ifa+l)  + S£  Ne  Nq  + £9 A/e 


(204) 


where  N^,N^  are  respectively  the  maximum  number  of  gores  and  the  maximum 
number  of  elements  per  gore  to  be  specified  in  any  case  included  in  the 
run.  Values  of  N^  up  to  and  including  10  elements  per  gore  are  permitted. 
The  dimension,  as  determined  from  Eq.  204  must  also  be  defined  in  the  DATA 
statement : 

DATA  LENGTH, NLY , NLYPl/xxx ,y , z/ 


where  also: 


y = total  number  of  distinct  layers  in  the  laminated  plate 
z = y+1 
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The  auxiliary  arrays  and  vectors  must  then  be  dimensioned  as  follows  to 
agree  with  the  DATA  statement: 


DIMENSION  CM  (y, 3, 3,10) , CMC(y,3,3),  H(y),  XEL(y,6),  XR(y),  Z(z) 


Second,  a job  control  instruction  external  to  the  FORTRAN  code  must  be 

provided  to  establish  a temporary  sequential-access  dataset  on  system  disk 

or  tape  storage.  The  dataset  must  be  allocated  418  single-precision  words 

per  record  and  a number  of  records  at  least  equal  to  N , the  maximum  number 

9 

of  gores  to  be  used  in  any  case  in  the  run.  The  dataset  must  be  assigned 
FORTRAN  unit  number*  20.  Job  control  instructions  must  also  be  provided  to 
access  the  pre-compiled  codes,  or  the  source  decks  must  be  provided  for  the 
general  software  modules  shown  in  Fig.  32. 

Third,  the  user  must  prepare  input  data  cards  in  accordance  with  the 
conventions  given  in  Table  10.  Note  that  cases  with  different  amounts  of 
mesn  refinement  (NGOR  N^ , NEPG  N^)  may  be  stacked  together,  and  that 
various  combinations  of  loading  may  be  applied  to  the  models.  Also,  note 
that  the  correct  ply  angle  9 for  a layer  of  the  laminate  is  opposite  in 
sign  to  the  wrapping  angle  (j>  which  would  be  specified  for  a filament-winding 
operation  to  fabricate  the  cone.  The  sign  difference  results  from  the 
irrangement  of  coordinate  systems,  as  shown  in  Fig.  33.  The  figure  also 
shows  how  the  stress  solution  printout  {a  O O a a } in  each  local 


y xy  xz  yz 


gore  plane  may  be  interpreted  approximately  as  shell  stresses  {a.  aQ  a. 

Oa  ^ 0^},  where represents  the  shell  surface  coordinate  system. 

Subroutine  AMMRC4  (Appendix  C)  is  a modification  of  subroutine  AMMRC3 
in  which  two  options  have  been  added  to  the  code.  First,  a fourth  mode  of 
applied  loading  is  permitted,  in  the  form  of  a total  side-force  S which  is 
distributed  over  that  portion  of  the  shell  surface  in  the  region  X<0,  i.e. 
n/2  < 0 < 3ir/2  (see  Fig.  31).  The  distribution  is  assumed  to  be  a pressure 
loading  normal  to  the  shell  surface,  uniform  with  respect  to  Z and  varying 
sinusoidally  with  0: 


(R)  t'j  2 .)  - -jo0  cos  P 


(205) 


See  footnote  on  page  91  for  definition  of  FORTRAN  unit  number. 
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where 


To  = 


2 s -/iy*i)z  ^ 

CKj  + R,)  Ll 


(206) 


The  side-foirce  is  intended  to  simulate  local  aerodynamic  or  blast  loading 
on  the  shell. 

Second,  a ring  of  stiffeners  may  be  added  to  the  free  tip  of  the  shell 
(Z=L,  R=R^)  to  simulate  the  presence  of  a metal  grip  fixture  such  as  would 
be  used  in  experimental  apparatus  to  introduce  the  tip  loads  M,F,V  into  the 
shell  structure.  Subroutine  STIF2,  a standard  assumed-displacement  engineer- 
ing beam  theory  element,  is  used  for  this  purpose.  Element  STIF2  incorporates 
cubic  bending  behavior  in  two  planes,  linear  axial  stretching,  and  linear 
axial  twisting  in  its  displacement  field.  Figure  34  illustrates  the  element, 
together  with  definitions  of  the  section  properties  which  are  required  to 
describe  it.  In  the  present  case,  the  element  cross  section  is  oriented  as 
shown  in  Fig.  35  for  the  purpose  of  calculating  inputs  for  the  cross  section 

inertias  I , 1 and  I . The  element  nodal  coordinates,  together  with 
yy  zz  yz 

addition  of  a sixth  degree  of  freedom  at  each  node  coupled  to  a stiffener, 
are  programmed  internally  in  the  AMMRC4  code. 

Subroutine  AMMRC4  is  prepared  for  execution  in  the  same  manner  as 
subroutine  AMMRC3 , with  the  following  exceptions.  First,  a slightly  larger 
dimension  for  the  FEABL-2  data  vector  may  be  needed  if  the  stiffener  option 
is  exercised.  Second,  the  conventions  for  input  data  cards  have  been  modified 
slightly,  as  shown  in  Table  11. 

Only  limited  verification  tests  of  subroutines  AMMRC3  and  AMMRC4  were 
conducted  because  the  performance  of  all  basic  modules  used  by  these  codes 
had  previously  been  verified  (see  Section  7).  Hence,  only  a brief  test 
problem  of  an  isotropic  cylinder  was  run  simply  to  verify  the  mesh-  and 
load-generation  algorithms.  After  verification,  a few  runs  of  laminated 
conical  shells  were  made  for  demonstration  purposes . Subroutines  AMMRC3 
and  AMMRC4  are  intended  primarily  for  future  verifications  involving  the 
comparison  of  computed  stresses  and  strains  with  experimental  data  currently 
being  sought  by  AMMRC  under  another  program. 
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Subroutines  AMMRC2,  AMMRC3  and  AMMRC4  have  been  programmed  to  be 
compatible  with  the  existing  FEABL-2  software  [6] . However,  these  codes 
may  be  modified  easily  for  compatibility  with  FEABL-5.  The  following 
actions  are  required  to  modify  each  subroutine.  First  labelled  COMMON  area 
SIZE  must  be  expanded  to: 


COMMON  /SIZE/  NET , NDT , NUT , NSP , I ODYN 


Second,  the  statement  IODYN=0  must  be  placed  near  the  beginning  of  the 
subroutine.  The  other  added  parameters  NUT, NSP  need  not  be  defined  for  the 
static  analyses  executed  by  the  AMMRC  subroutines.  Finally,  all  instructions 
which  generate  prescribed  displacements  and  prescribed  nodal  forces  must  be 
moved  from  their  present  locations  into  positions  between  the  calls  to  sub- 
routines SETUP  and  ORK,  and  this  must  be  done  without  changing  the  order  in 
which  these  instructions  occur  relative  to  each  other. 


8.5  Subroutine  VIBEPT 


Subroutine  VIBEPT  was  programmed  to  verify  integration  of  the  FEABL-5 
dynamic  analysis  software.  The  code  generates  and  computes  eigenvalues  for 
a quarter-model  of  a simply  supported  rectangular  plate.  After  the  eigen- 
value analysis  has  been  completed,  a transient  response  analysis  problem  is 
executed.  Subroutine  VIBEPT  was  used  to  perform  the  test  analyses  described 
in  Subsections  7.5.2  and  7.5.3.  The  general  flow  chart  given  in  Fig.  36 
illustrates  the  complete  sequence  of  operations  required  in  a FEABL-5  dynamic 
analysis.  Subroutine  VIBEPT  is  listed  in  Appendix  D. 

The  following  user  actions  are  required  to  execute  the  analysis.  First, 
vector  and  array  dimensions  must  be  properly  defined  in  the  dummy  main  program 
Suggested  dimensions  for  the  FEABL-5  data  vector  range  from  2,500  words  for 
a 3x3-element  mesh,  to  52,800  words  for  a 9x9-element  mesh,  as  detailed  in 
Table  12.  Whatever  dimension  is  chosen  must  also  be  defined  in  the  DATA 


DATA  NLY , NSPACE , LENGTH/y , P , xxx/ 


y = total  number  of  distinct  layers  in  the  laminated  plate. 

P = total  number  of  degrees  of  freedom  to  be  included  in  the 
subspace  for  eigenvalue  iteration  (see  Subsection  7.2). 
xxx  = length  of  the  data  vector. 


Also,  the  auxiliary  arrays  and  vectors  must  be  dimensioned  as  follows  to 
agree  with  the  DATA  statement: 


DIMENSION  CMC (y, 3, 3) , H (y) , DENS (y) , XR(y),  XEL(y,6),  Z(z) 

DIMENSION  ELK(q),  EMS (q) , DELK(q),  DEMS (q) 

2 

DIMENSION  EV(P  ),  T(P,P),  AM(P,P),  ICONV(P),  ASQ(P) 


where  z=y+l  and  q=P(P+l)/2.  Second,  the  user  must  prepare  input  data  cards 
in  accordance  with  the  conventions  given  in  Table  13.  Subroutine  VIBEPT 
does  not  require  any  temporary  datasets.  However,  job  control  instructions 
must  be  provided  to  access  the  pre-compiled  codes,  or  the  source  decks  must 
be  provided  for  the  general  software  modules  shown  in  Fig.  37. 

Results  of  the  transient  response  verification  which  were  run  with 
subroutine  VIBEPT  were  presented  in  Subsection  7.5.3.  Analyses  of  a 3x3- 
element  model  were  conducted  for  two  cases  of  sinusoidally  varying  transverse 
load,  one  case  in  which  the  plate  was  given  initial  displacement,  and  one 
case  in  which  the  plate  was  given  an  initial  velocity.  These  various  condi- 
tions are  established  by  subroutine  INPUT,  which  was  programmed  as  a temporary 
code  to  establish  load  time-history  and  initial  condition  data  for  FEABL-5 


subroutines  TSTEP  and  ELAPSE.  These  subroutines  are  listed  in  Appendix  E. 
The  version  of  subroutine  INPUT  given  in  the  appendix  establishes  zero 
initial  conditions  and  the  sinusoidal  load  time-history  with  a spatial 
distribution  corresponding  to  the  first  natural  mode  of  the  plate.  Sub- 
routine TSTEP  will  be  modified  to  accept  general  load  time-histories  prior 
to  its  documentation  in  the  FEABL-5  user's  guide. 


SECTION  9 


DISCUSSION  AND  CONCLUSIONS 


9.1  Summary 

This  report  has  presented  the  results  of  an  extension  of  previous 
investigations  into  the  formulation  of  assumed-stress  hybrid  finite  elements 
for  bending  of  multilayer  laminated  plates  and  shells.  The  previous  investiga 
tions  resulted  in  a quadrilateral  thick-plate  element  capable  of  reproducing 
the  severe  cross-section  warping  which  results  when  an  extremely  thick  plate 
consisting  of  only  a few  elastically  dissimilar  layers  is  subjected  to  bend- 
ing. During  the  current  investigation,  a triangular  version  of  this  element 
was  formulated,  and  attempts  were  made  to  modify  the  existing  quadrilateral 
element  by  incorporating  an  assumed  stress  distribution  which  would  exactly 
satisfy  traction-free  boundary  conditions  along  one  of  the  element's  lateral 
edges.  The  candidate  traction-free  stress  distributions  which  were  investi- 
gated in  this  study  were  unsuccessful  in  that  they  led  to  internal  kinematic 
instabilities  in  the  element  stiffness  matrix.  Also,  study  of  the  practical 
aspects  of  this  element  family  revealed  that  its  employment  for  analysis 
of  typical  advanced  fiber  composite  laminates  (which  generally  consist  of 
ten  or  more  elastically  distinct  layers)  would  require  more  core  memory 
than  is  generally  available,  even  in  today's  large  digital  computer  systems. 

A second  family  of  multilayer  plate  bending  elements  was  subsequently 
investigated  in  order  to  provide  similar  analysis  capabilities  without  the 
storage  penalty.  The  new  family  is  based  on  stress  distributions  derived 
from  assumed  in-plane  strains,  and  thus  does  not  reproduce  severe  cross 
section  warping.  However,  the  transverse  shear  effects  associated  with  the 
behavior  of  moderately- thick  plates  are  included,  and  these  elements  are 
thus  applicable  to  both  thin  and  moderately-thick  plates.  The  software 
which  resulted  from  this  part  of  the  study  has  been  programmed  in  the  form 
of  independent  modules:  subroutines  MLP3K(Q/T)  and  MLTPK  which  compute 

the  stiffness  matrices  for,  respectively,  one  quadrilateral  and  two  tri- 
angular elements,  subroutine  MLP3S  which  computes  element  stresses  and 
strains  from  given  nodal  displacements,  and  subroutine  MLP3M  which  computes 
both  lumped  and  hybrid-rational  element  mass  matrices  for  dynamic  analyses. 
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The  modularity  of  these  subroutines  permits  them  to  be  adapted  to  any 
finite-element  code  with  little  or  no  reprogramming  required.  In  the  present 
investigation,  the  subroutines  were  combined  with  the  existing  ASRL  FEABL-2 
code,  which  contains  the  global  algorithms  for  assembly  and  solution  of 
finite-element  equation  systems  and  which  is  also  modular.  In  the  course 
of  the  investigation,  FEABL-2  was  modified  for  compatibility  for  several 
additional  dynamic  analysis  modules  which  were  adapted  from  concepts  presented 
in  the  previous  investigations.  The  modified  software,  designated  FEABL-5, 
will  be  documented  in  detail  in  a separate  user's  guide.  The  dynamic  analysis 
software  consists  of  two  major  modules.  The  first  executes  an  eigenvalue 
analysis  to  compute  the  natural  frequencies  and  mode  shapes  associated  with 
an  assembled  finite-element  model  of  a structure.  This  module  is  based  on 
the  Subspace  Iteration  Method,  which  has  proved  to  be  an  efficient  approach 
to  modal  analysis  where  the  primary  interest  is  to  obtain  the  first  few 
modes  (say  up  to  20  percent  of  the  total  number  represented  in  the  model) . 

The  second  module  executes  a transient  response  analysis  based  on  the  Modal 
Superposition  Method  (MSM) . Arbitrary  damping  factors  (up  to  critical)  for 
each  mode,  arbitrary  initial  displacement  and  velocity  conditions,  and 
arbitrary  applied  force  time-histories  may  be  input.  The  MSM  was  chosen 
rather  than  a direct  time-integration  finite-difference  approach  for  the 
following  reasons: 

1.  The  natural  frequencies  and  mode  shapes  required  for  the  MSM  are 
likely  to  be  of  interest  anyway  for  various  engineering  purposes. 

2.  Because  an  analytic  solution  is  used  in  the  MSM,  one  can  elect  to 
calculate  response  information  at  any  time  instant  or  sequence  of 
time  instants  apart,  spaced  evenly  or  unevenly.  There  is,  in 
principle,  no  restriction  on  the  size  of  the  time  interval  between 
"sampling"  instants,  except  for  that  cited  next  in  item  3.  Further, 
this  type  of  solution  does  not  introduce  frequency  distortion  or 
false  damping;  such  detrimental  effects  are  often  encountered  in 
direct  timewise  finite-difference  solution  methods. 

3.  In  the  MSM,  the  time-step  size  or  sampling  interval  can  easily  be 
varied  arbitrarily  from  one  step  to  the  next.  This  is  particularly 
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convenient  for  modelling  load  time-histories  in  which  periods  of 
slowly  and  rapidly  changing  loads  are  interspersed.  Of  course, 
one  must  insure  that  an  adequate  number  of  modes  has  been  included. 

It  should  be  realized,  however,  that  if  the  spatial  distribution  of  the 
externally-applied  loading  varies  significantly  with  time,  one  is  faced 
with  a considerable  amount  of  calculation  in  order  to  evaluate  the  generalized 
force  acting  on  each  mode;  this  type  of  added  computational  burden  is  not 
nearly  as  severe  when  one  employs  the  generalized-displacement  direct  time- 
wise  solution  procedure. 

A further  modification  of  the  new  family  of  plate  elements  was  also 
studied.  In  this  case,  an  extra  partial  layer  was  added  to  the  laminate  to 
simulate  the  presence  of  an  integral  stiffener  offset  from  the  plate  neutral 
axis.  This  formulation  was  compared  with  the  conventional  approach  to 
stiffened  plates  by  means  of  another  existing  module,  subroutine  STIF2  which 
is  a separate  assumed-displacement  model  of  a stiffener. 

Finally,  the  various  modules  described  above  were  put  together  in 
several  combinations  to  form  applications  programs  for  different  purposes. 
These  applications  programs  included  analyses  of  a rectangular  plate  with 
circular  hole,  subjected  to  four-point  bending;  a conical  shell  subjected 
to  tip  shear,  axial  load,  and  bending  and  a distributed  side-force;  and  an 
eigenvalue  and  transient  response  analysis  of  a simply-supported  flat  plate. 

9.2  Discussion  of  Results 

Comparative  performance  tests  of  the  thick  and  moderately-thick  plate 
elements  were  conducted  to  assess  the  abilities  of  these  elements  to  reproduce 
the  bending-stretching  behavior  which  is  expected  to  occur  in  multilayer 
laminated  plates  and  shells.  The  assessments  were  made  by  comparison  with 
independent  analytical  solutions  obtained  by  other  investigators.  The  thick- 
plate  elements  were  found  to  represent  accurately  the  severe  cross  sectional 
warping  which  occurs  in  short-span  laminates  (span/thickness  ratio  S/t  < 10) 
which  have  less  than  4 distinct  layers.  Accurate  results  were  obtained  for 
very  thick  plates  (S/t  =4).  On  the  other  hand,  these  elements  exhibited 
the  distinct  disadvantage  that  the  requirements  for  total  degrees  of  freedom, 
core  memory  and  computation  time  rose  extremely  rapidly  as  the  number  of 
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layers  in  the  laminate  increased.  Four  layers  appeared  to  constitute  a 
practical  upper  limit,  given  the  capabilities  of  current  large  digital 
computing  facilities  (e.g.,  IBM  S-370/168  with  125,000  decimal  words  of 
core  memory  or  CDC-6600  with  300,000  octal  words). 

The  moderately- thick  plate  elements  were  observed  to  be  much  more 
efficient  in  this  respect.  For  these  elements,  the  total  number  of  degrees 
of  freedom  in  a model  is  independent  of  the  number  of  layers  in  the  laminate, 
as  is  the  core  memory  requirement,  while  the  computation  time  was  observed 
to  increase  insignificantly  as  the  number  of  layers  was  increased.  Demonstra- 
tion analyses  of  conical  shells  with  several  hundred  degrees  of  freedom  and 
up  to  25  distinct  layers  required  only  on  the  order  of  1 CPU  minute  for  mesh 
generation  and  a complete  static  stress/strain  solution.  The  technical 
performance  of  these  elements  was  observed  to  be  poor  for  thick  plates  (S/t=4) 
with  few  layers,  as  expected.  However,  the  results  for  thick,  many-layered 
plates  (S/t=4,  7 layers)  and  the  results  for  moderately  thick  to  thin  plates 
(S/t>10)  were  observed  to  be  accurate.  (In  these  cases,  the  exact  stress 
solutions  tend  strongly  to  approach  lamination  theory,  and  the  cross- 
sectional  warping  effects  become  much  less  significant.)  It  is  important 
to  recall  that  advanced  composite  plates  and  shells  generally  consist  of  ten 
or  more  distinct  layers  and  have  S/t>10. 

The  attempt  to  formulate  a special  traction-f ree-edge  (TFE)  thick-plate 
element  was  motivated  by  a concern  about  solution  accuracy  in  the  vicinity 
of  free  lateral  edges  of  plates  and  shells.  Past  experience  with  assumed- 
stress  hybrid  elements  has  shown  that  the  TFE  modification  is  sometimes 
(but  not  always)  required  to  obtain  accurate  stresses  near  free  edges 
associated  with  geometrical  details  (e.g.,  a circular  hole  in  a rectangular 
plate) . As  has  been  mentioned,  the  attempt  to  formulate  a TFE  thick-plate 
element  proved  unsuccessful.  A TFE  moderately- thick  plate  element  was  also 
considered  briefly.  However,  an  investigation  of  the  formulation  of  this 
element  quickly  revealed  that  such  a modification  was  mathematically  impossible. 
The  applications  program  mentioned  previously  (analysis  of  a plate  with  circular 
hole  in  four-point  bending)  was  put  together  specifically  to  assess  the  poten- 
tial inaccuracies  which  might  occur  when  the  non-specialized  hybrid  elements 
were  used  to  model  geometrical  details  of  this  type.  Fortunately,  the 


computed  solutions  were  observed  to  reproduce  all  free-edge  conditions  with 
good  accuracy,  for  both  the  thick  and  moderately- thick  plate  elements.  The 
computed  solutions  were  also  observed  to  reproduce  accurately  the  bending 
stress  concentration  near  the  hole,  as  was  shown  by  comparison  with  inde- 
pendently obtained  experimental  data. 

Attention  was  then  focussed  strictly  on  the  moderately-thick  plate 
elements,  since  these  appeared  to  be  the  more  practical  analysis  tools  for 
composite  laminates.  Parametric  analyses  of  unbalanced  two-layer  (+0) 
simply  supported  rectangular  plates  were  conducted  and  compared  with  an 
independent  analytical  solution  to  obtain  additional  performance  data.  The 
quadrilateral  element,  MLP3K(Q),  was  observed  to  produce  results  within  3 
percent  of  the  exact  solutions  for  stresses  and  displacements  with  the  ply 
angle  varied  parametrically,  15°  9 _<  75°.  On  the  other  hand,  the  triangle 

element  MLP3K (T)  proved  to  be  too  stiff.  Errors  of  30  percent  were  observed 

o 

when  this  element  was  subjected  to  a convergence  tests  for  the  case  of  a +45 
square  plate.  These  results  led  to  the  creation  of  element  MLTPK,  a triangle 
element  similar  to  MLP3K(T),  but  with  fewer  assumed-stress  parameters  to  give 
the  element  additional  flexibility.  The  convergence  study  with  element  MLTPK 
resulted  in  errors  of  10  percent,  i.e.  within  the  bounds  of  reasonable 
engineering  accuracy. 

In  another  series  of  tests  of  element  MLP3K(Q),  the  ply  angle  was  fixed 
at  9=45°,  while  the  plate  shape  and  number  of  elements  in  each  direction 
were  varied  to  produce  different  aspect  ratios  (length/width) . The  results 
of  these  tests  indicated  that  no  apparent  inaccuracies  resulted  from  -.^Dect 
ratios  as  large  as  10.  These  tests  were  conducted  with  the  critical  element 
stiffness  computations  carried  out  in  double-precision  arithmetric  (approxi- 
mately 15  significant  decimal  figures).  However,  severe  aspect  ratio  degrada- 
tions were  observed  in  MLP3K(£>)  elements  with  aspect  ratios  of  7 when  these 
calculations  were  carried  out  in  single-precision  (approximately  7 decimal 
significant  figures) . The  latter  results  were  obtained  during  some  of  the 
trail  analyses  of  the  plate  with  a circular  hole. 

In  another  series  of  tests,  element  MLP3K(Q)  was  modified  to  incorporate 
an  extra  partial  layer  simulating  offset  integral  stiffeners.  The  modified 
MLP3K(Q)  and  the  unmodified  MLP3K(Q)  combined  with  separate  assumed-displacement 


stiffener  elements  were  used  to  analyze  a simply  supported  square  plate 
with  two  integral  stiffeners  parallel  to  each  pair  of  edges.  A Fourier 
analysis  was  also  carried  out  to  provide  an  independent  analytical  solution 
for  comparison.  The  results  of  these  tests  showed  that  the  combination  of 
separate  stiffeners  with  unmodified  MLP3K(Q)  elements  gave  accurate  answers, 
while  the  modified  MLP3K(Q)  gave  poor  answers.  The  latter  result  was, 
tentatively,  attributed  to  the  violation  of  traction- free  conditions  on  the 
lower  surface  of  the  plate,  a situation  caused  by  the  modified  stress  assump- 
tions . 

Tests  of  the  dynamic  analysis  modules  were  conducted  primarily  for 
verification  and  demonstration  of  software  compatibility.  The  conclusions 
of  previous  investigators  about  the  efficiency  of  the  Subspace  Iteration 
Method  of  eigenvalue  analysis  were  reconfirmed,  with  all  tests  limited  to 
the  simplest  of  three  previously  proposed  methods  for  generating  initial 
estimates  of  the  natural  mode  shapes.  However,  three  additional  interesting 
observations  resulted  from  these  tests.  First,  the  converged  eigenvalues 
were  found  to  be  quite  sensitive  to  modelling  detail.  Errors  which  could 
not  be  reduced  by  continued  iteration  were  observed  for  eigenvalues  correspond- 
ing to  mode  shapes  with  spatial  distributions  having  variations  more  rapid 
than  the  capability  of  the  MLP3K(Q)  element  interpolation  functions.  Second, 
the  rates  of  convergence  of  the  eigenvalues  were  found  to  have  some  sensitivity 
to  the  number  of  modes  included  in  the  subspace.  However,  the  asymptotic 
convergence  rate  theory  developed  by  other  investigators  was  found  to  give 
poor  predictions  for  required  computation  times.  Third,  comparisons  of 
analyses  performed  with  lumped  and  hybrid-rational  mass  matrices  showed  that 
the  hybrid- rational  approach  gave  better  results  for  the  first  mode,  while 
the  lumped-mass  approach  gave  better  results  for  the  second  and  higher  modes. 
Although  the  errors  in  both  cases  were  well  within  reasonable  engineering 
accuracy,  this  result  is  a surprising  reversal  of  the  behavior  commonly  found 
when  lumped  and  consistent  mass  models  are  compared  in  analyses  based  upon 
assumed-displacement  elements. 

The  Modal  Superposition  Method  software  module  was  verified  in  a brief 
series  of  tests  in  which  a simply  supported  square  plate  was  modelled. 

Accurate  results  were  obtained  in  four  test  problems:  free  vibration  of  the 
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plate  with  a prescribed  initial  displacement  field,  free  vibration  with  a 
prescribed  initial  velocity  field,  and  two  cases  of  steady-state  sinusoidal 
forced  vibration. 

Finally,  an  informal  assessment  of  the  effectiveness  of  the  modular 
concept  was  conducted  during  the  applications  program  phase  of  the  study. 

This  phase  of  the  work  involved  several  analysis  modifications  which  were 
requested  by  AMMRC , and  which  were  performed  with  a minimum  amount  of 
reprogramming  by  ASRL  personnel.  More  significantly,  other  similar  modifica- 
tions were  carried  out  by  AMMRC  personnel  in  parallel  with  learning  to  use 
the  various  modules.  These  latter  modifications  were  generally  accomplished 
rapidly,  after  receipt  of  brief  advice  and  suggestions  from  ASRL.  One  such 
modification  has  resulted  in  the  coupling  of  the  FEABL  software  with  an 
AMMRC  in-house  program  for  automatic  generation  of  finite-element  grids 
for  structures  with  boundaries  of  arbitrary  shape. 

9.3  Conclusions  and  Recommendations 

Several  useful  conclusions  can  be  drawn  from  the  results  discussed 
above,  as  follows: 

1.  The  family  of  moderately- thick  plate  elements  appears  to  suit  the 
analysis  needs  for  practical  advanced  fiber  composite  plates  and 
shells  better  than  does  the  family  of  thick-plate  elements.  It 

is  recommended  that  element  MLP3K(Q)  be  retained  in  a primary  role, 
with  element  MLTPK  serving  a secondary  role  as  a mesh-expander. 

2.  A traction-f ree-edge  special-purpose  element  does  not  appear  to  be 
required  for  plates  with  circular  cutouts;  numerical  accuracy  of 
the  unmodified  elements  has  been  demonstrated  in  tests  with  free- 
edge  geometrical  details. 

3.  The  family  of  moderately-thick  plate  elements  appears  to  have 
some  performance  limitations  in  terms  of  aspect  ratio  (length/ 
width) . Element  aspect  ratios  of  less  than  5 are  recommended  for 
the  elements  in  present  (single-precision)  form,  while  aspect 
ratios  of  7 to  10  appear  to  be  acceptable  if  the  element  stiffness 
subroutines  are  converted  to  double-precision  arithmetic. 

4.  The  current  procedure  of  modelling  integrally  stiffened  plates  by 
combining  separate  plate  and  stiffener  elements  is  acceptable  for 
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engineering  results.  Modification  of  the  plate  elements  to 
incorporate  integral  stiffeners  proved  to  be  unsatisfactory. 
However,  all  of  the  possible  approaches  to  the  latter  type  of 
formulation  have  not  yet  been  exhausted. 

The  combination  of  the  Subspace  Iteration  Method  and  the  Modal 
Superposition  Method  appears  to  be  the  most  efficient  approach 
for  the  type  of  dynamic  analysis  required  to  evaluate  the  small 
displacement  linear-elastic  transient  responses  of  advanced 
composite  plates  and  shells. 

6.  Modularization  of  the  finite-element  software  has  proved  to  be  of 
great  benefit  in  terms  of  both  ease  of  transfer  from  the  develop- 
ing to  the  using  organization  and  flexibility  in  modification  for 
future  analysis  tasks. 

In  addition  to  the  above  conclusions,  several  recommendations  for 
possible  future  developments  are  presented.  First,  additional  carefully- 
controlled  tests  of  combinations  of  elements  MLP3K(Q)  and  MI.TPK  should  be 
carried  out.  Element  MLTPK  has  been  identified  as  a good  "mesh-expander", 
i.e.  an  element  which  can  be  used  to  transist  between  coarse  grids  of 
quadrilaterals  in  regions  with  low  stress  gradients  to  fine  grids  in 
regions  with  high  stress  gradients.  The  numerical  errors  in  such  analyses 
will  probably  be  lower  than  the  10  percent  level  observed  in  tests  of  models 
consisting  only  of  MLTPK  elements.  However,  this  is  merely  a conjecture 
until  confirmed  by  numerical  experiment. 

Second,  it  appears  to  be  worthwhile  to  investigate  the  possibility 
of  modifying  element  MLP3K(Q)  to  incorporate  the  laminate  edge  effects  which 
are  known  to  occur  near  free  lateral  edges.  The  edge  effects  spoken  of  here 
do  not  refer  to  the  TFE-type  modifications  discussed  previously.  Rather, 
the  intent  is  to  model  properly  the  peaks  in  the  interlaminar  peel  and 
shear  stresses  which  appear  near  the  lateral  edges  of  laminates.  These 
stresses  constitute  a departure  from  simple  lamination  theory,  and  the 
presence  of  the  peel  stress  requires  displacement  interpolations  which 
permit  the  transverse  displacement  to  vary  through  the  laminate  thickness. 
Hence,  extra  degrees  of  freedom  are  required  along  one  edge  of  the  modified 
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element,  as  well  as  revised  interpolations  for  the  assumed  stress  field. 

This  possibility  is  judged  to  be  worth  investigation  because  it  offers 
the  potential  of  a natural  incorporation  of  the  important  edge  effects 
within  the  finite- element  method. 

Third,  it  is  recommended  that  additional  alternatives  for  integrally 
stiffened  plate  elements  be  investigated.  Although  the  present  method  of 
combining  assumed-displacement  stiffeners  with  unstiffened  plate  elements 
is  workable,  it  introduces  the  practical  inconvenience  that  finite-element 
grids  must  be  arranged  so  that  all  stiffeners  lie  along  edges  of  plate 
elements.  This  can  result  in  an  over-refined  mesh  in  areas  of  low  stress 
gradient  where  many  stiffeners  are  placed.  The  principal  alternatives  to 
be  examined  are  other  modified  plate  stress  assumptions  and  rational  trans- 
formations of  separate  stiffener  elements  based  on  the  plate  element  dis- 
placement interpolation. 

Fourth,  the  development  of  other  special-purpose  elements  based  on  the 
assumed-stress  hybrid  method  may  be  useful  for  advanced  analyses  of  real 
composite  structures.  One  possible  example  is  an  element  which  includes 
within  its  interior  a fastener  detail  or  cutout.  An  element  of  this  type 
has  already  been  developed  and  applied  successfully  to  problems  involving 
plane  stress  analysis  and  fracture  mechanics  analysis  of  single- layer 
isotropic  media. 

Finally,  it  must  be  recognized  that  even  the  Subspace  Iteration  Method 
will  tax  the  capacity  of  current  digital  computers  for  eigenvalue  analyses 
of  the  very  detailed  finite-element  models  which  may  be  required  to  assess 
the  dynamic  behavior  of  real  composite  structures.  Therefore,  it  is 
recommended  that  application  of  the  Component  Mode  Synthesis  Method  [27]  to 
these  analyses  be  investigated.  The  Component  Mode  Synthesis  Method  is 
a rational  approach  to  substructuring  a finite-element  model  for  dynamic 
analysis.  The  necessary  global  software  for  these  computations  is  currently 
being  developed  under  another  contract,  and  will  be  implemented  in  the  near 
future  as  an  additional  FEABL-5  module. 
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TABLE  3 

A COMPARISON  OF  ANALYSES  FOR  THE  STIFFENED 
SIMPLY  SUPPORTED  PLATE 


CENTER  LATERAL  DEFLECTION 


t 

s 

U) 

s 

MESH 

RITZ 

METHOD 

INTEGRALLY 

STIFFENED 

NON INTEGRALLY 
STIFFENED 

2x2 

0.  300 

0.300 

4x4 

0.297 

0.297 

0.0 

0.0 

0. 

296 

6x6 

0.297 

0.297 

8x8 

0.297 

0.  297 

2x2 

0.299 

0.300 

4x4 

0.297 

0.297 

0.01 

0.03 

0. 

296 

6x6 

0.297 

0.297 

8x8 

0.297 

0.297 

2x2 

0.298 

0.299 

4x4 

0.295 

0.297 

0.03 

0.03 

0. 

295 

6x6 

0.295 

0.297 

8x8 

0.294 

0.297 

0.1 



0.03 

2x2 

4x4 

6x6 

8x8 

0. 

293 

0.  291 
0.284 

0.  270 
0.251 

0.297 

0.295 

0.295 

0.295 

TABLE  4 

EFFECT  OF  ARITHMETIC  PRECISION  ON  THE 
JACOBI  ITERATION  METHOD 


ANALYTICAL  SOLUTION  (1) 

— 

JIM  SOLUTION  IN 
SINGLE  PRECISION  (2) 

JIM  SOLUTION  IN 
DOUBLE  PRECISION  (2) 

VALUE 

% ERROR 

VALUE 

% ERROR 

CANTILEVER: 

U)  = (1.875)  2 

/ 4 

EI/mL 

249.4 

249.4 

0 

249.4 

0 

W2 

(4.694) 2 

y ~4 

EI/mL 

1,  563 

1,  564 

.064 

1,564 

.064 

W3 

(7.855) 2 

J " 4 

EI/mL 

4,  376 

4,392 

.366 

4,  392 

.366 

II 

3 

(11. oo)2 

/ 4 

EI/mL 

8,576 

8,676 

1.  17 

8,676 

1.17 

W5 

(14. 14)2 

/ 7 4 

EI/mL 

14, 180 

14,400 

3.12 

14,400 

3.12 

W6  = 

(17. 28)2 

/ " 4 

EI/mL 

21,180 

23,920 

12.96 

23,920 

12.96 

W7  = 

(20.42) 2 

/ 7 ' 4 

EI/mL 

29,580 

34,990 

18.32 

34,990 

18.32 

II 

00 

3 

(23. 56)2 

/ 4 

EI/mL 

39, 380 

50,740 

28.90 

50, 740 

28.90 

8 

II 

(26.70) 2 

/ 7 4 

EI/mL 

50,580 

72,080 

42.50 

72,080 

42.50 

U)  = 

10 

(29.85) 2 

/ 4 

EI/mL 

63,180 

106,000 

67.70 

106,000 

67.70 

UNRESTRAINED: 

0)^  = 0 

0 

32.46 



0.00184 

— 

^2  = 

0 

4 

EI/mL 

0 

-56.02 

— 

0.0348 

— 

u»3  = 

(1 . 506tt)  2 

1,587 

1,587 

.063 

1,588 

.063 

II 

3 

(2.5007T)2 

/ . 4 

EI/mL 

4,  375 

4,388 

.298 

4,388 

.298 

I 

I j 

I 

I 


J. 


Notes:  (1)  Analytical  solutions  from  Ref.  20.  Numerical  values  are  for 

5 2 

the  case  /EI/in  = 1.134901x10  in  /sec  and  L = 40  inches. 

Solutions  are  in  units  of  rad/sec. 

(2)  JIM  solutions  obtained  from  12-DOF  finite-element  model 

-4 

(5  elements)  and  with  tolerance  c = 10  . Arithmetic  precisions 

for  IBM  S- 370/168  are:  Single  = 32  BITS  =7.2  decimal  significant 

figures,  Double  = 64  BITS  = 15  decimal  significant  figures. 
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TABLE  5 


EFFECT  OF  MODELLING  DETAIL  ON  THE  JACOBI 
ITERATION  METHOD 


ANALYTICAL  SOLUTION  (1) 


JIM  SOLUTION  (2) 
12- DOF  MODEL 

JIM  SOLUTION  (2) 
50-DOF  MODEL 

VALUE 

% ERROR 

VALUE 

% ERROR 

249.4 

0 

249.4 

0 

1 , 564 

.064 

1,563 

0 

4,392 

. 366 

4,  376 

0 

8,676 

1.17 

8,576 

0 

14, 400 

3.12 

14,180 

0 

23,920 

12.96 

21,180 

0 

Notes:  (1)  Analytical  solutions  from  Ref.  20.  Numerical  values  are  for 

. 5 2 

the  case  /EI/m  = 1.134901x10  in  /sec  and  L =40  inches. 
Solutions  are  in  units  of  rad/sec. 

(2)  JIM  solutions  in  double-precision  arithmetic  (64  BITS  = 15 
decimal  significant  figures)  on  IBM  S-370/168. 


♦MODELED  QUARTER  PLATE  INCLUDING  IN-PLANE  DEGREES  OF  FREEDOM 


TABLE  7 


THEORETICAL  EFFECT  OF  SUBSPACE  SIZE  ON 


EIGENVALUE  CONVERGENCE  RATES 


0.90381 


Analytical  solutions  from  Ref.  25 


SIZE 

OF  SUBSPACE,  P = 

6 

8 

9 

10 

12 

.00346 

.00160 

.00160 

.00160 

.00119 

.0865 

.0400 

.0400 

.0400 

.0297 

.0865 

.0400 

.0400 

.0400 

.0297 

.280 

.130 

.130 

.130 

.0963 

.585 

.270 

.270 

.270 

.201 

.585 

.270 

.270 

.270 

.201 

CONVERGENCE  RATES: 
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INPUT  CONVENTIONS  FOR  SUBROUTINE 


AEROELASTIC  AND  STRUCTURES  RESEARCH  LABORATORY 

JOB  TITLE SUBROUTINE  AMM  R.C.  3 ENGINEER  d R&  JUf-eP. 

JOB  NO  826!*)  PROGRAMMER  O RRldyeR.  /FRE NCR  DAT  £ 12.  M *7  1976  SHEET  1 0 


INPUT  CONVENTIONS  FOR  SUBROUTINE  AMMRC3 
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TABLE  12 

VARIATION  OF  REQUIRED  CORE  STORAGE  WITH 
PROBLEM  SIZE  FOR  SUBROUTINE  VIBEPT 


MESH 

SIZE 

TOTAL 

DOF 

NUMBER  OF  DOF 
INCLUDED  IN 
SUBSPACE (P) 

LENGTH  OF 
FEABL-%  DATA 
VECTOR  (xxx ) 

REGION  SIZE 
SIZE* 

3x3 

80 

9 

2,500 

120K 

5x5 

180 

9 

9,900 

148K 

7x7 

320 

9 

24,700 

206K 

9x9 

500 

9 

49,300 

302K 

9x9 

500 

10 

50,500 

310K 

9x9 

500 

12 

52,800 

320K 

♦Region  size:  total  core  storage  required  for  problem  data  and  object 

code  on  IBM  S-370/168,  using  IBM  F0RTRAN-G1  and  FORTRAN- H (0)  compilers 

Region  sizes  are  given  in  KBYTES  (1  KBYTE  = 1,0241Q  single-precision  | 

words)  . 
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FIG.  1 POSITIVE  CONVENTION  FOR  SINGLE- ROTATION  AXIS  TRANSFORMATION 
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FIG.  2 PHYSICAL  SYMMETRIES  IN  A TYPICAL  FIBER-COMPOSITE  PLY 
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FIG.  23  SUMMARY  OF  FEABL-2  AND  FEABL-5  SOFTWARE  MODULES 
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Quadrant  Isolated  for  Finite-Element  Analysis,  with  Physical 
and  Symmetry  Boundary  Conditions  Indicated 


FOUR-POINT  BENDING  EXPERIMENT  ON  PLATE  WITH  CIRCULAR  HOLE 


FIG.  25  SCALE  PLANES  OF  TYPICAL  MESHES  GENERATED  BY  SUBROUTINE  AMMRC2 
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FIG.  35  ORIENTATION  OF  GRIP  RING  STIFFENERS  ASSUMED  IN  SUBROUTINE  AMMRC4 
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FIG.  37  MODULAR  ORGANIZATION  FOR  EXECUTION  OF  SUBROUTINE  VIBEPT 
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APPENDIX  B 


SUBROUTINE  AMMRC3 

This  appendix  contains  FORTRAN- IV  listings  of  applications  subroutine 
AMMRC3  and  the  dummy  MAIN  program  in  which  the  correct  dimensions  of 
variables  are  established. 


C MAIN  FOR  AMMRC3 

DIMENSION  RE(IOOOO),  IN(ICOOO) 

DIMENSION  CM(  1 ,3,  3,  10)  , CMC (1*3*3) f H(L),  XELll.ft),  XR(l),  Z < 1 J 
FCUIVALENCE  ( R1"  (1)  , IN(1)) 

DATA  I.  ENGTHfNLVfNLYPl/lOOOOtl  ,2/ 

CALL  AMMRC  )( LENGTH, NLY ,NLYP1 . RE  , IN  ,C“ .CMC .H.XEL.XR.Z) 

STOP 

END 
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APPENDIX  C 


SUBROUTINE  AMMRC4 

This  appendix  contains  FORTRAN- IV  listings  of  applications  subroutine 
AMMRC4  and  the  dummy  MAIN  program  in  which  the  correct  dimensions  of 
variables  are  established. 
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APPENDIX  D 
SUBROUTINE  VIBEPT 

This  appendix  contains  FORTRAN- IV  listings  of  applications  subroutine 
VIBEPT  and  the  dummy  MAIN  program  in  which  the  correct  dimensions  of 
variables  are  established. 
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Office  of  Secretary  of  Defense,  Office  of  the  Director  of  Defense  Research 
and  Engineering 

ATTN:  Mr.  J.  Persh,  Staff  Specialist  for  Materials  and  Structures 

The  Pentagon,  Washington,  D.C.  20301 

Commander,  U.S.  Army  Materiel  Command 

ATTN:  AMCRD-TT,  Dr.  R.  Zentner 

5001  Eisenhower  Avenue,  Alexandria,  VA  22333 

Ballistic  Missile  Defense  Program  Office,  ABMDA/W  (Provisional) 

ATTN:  DACS-BMT,  Mr.  C.  McLain 

DACS-BMT,  Mr.  V.  Kupelian 

Commonwealth  Bldg.,  Room  1100,  1300  Wilson  Blvd. , Arlington,  VA  22209 

Director,  Ballistic  Missile  Defense  Advanced  Technology  Center 
ATTN:  ATC-M,  Mr.  M.  Whitfield 

ATC-M,  Dr.  D.  Harmon 
ATC-X,  Mr.  W.  Davis 
P.0.  Box  1500,  Huntsville,  AL  35807 

Commander,  Ballistic  Missile  Defense  Systems  Command 
ATTN:  BMDSC-TEN,  Mr.  N.J.  Hurst 

P.O.  BOX  1500,  Huntsville,  AL  35807 

Director,  Defense  Nuclear  Agency 
ATTN:  SPAS,  Mr.  J.F.  Moulton,  Jr. 

SPAS,  Mr.  M.  Rubenstein 
Washington,  D.C.  20305 

Office  of  Chief  of  Research  Development  and  Acquisition, 

Department  of  the  Army 

ATTN:  DAMA-CSS,  Dr.  J.  Bryant 

Washington,  D.C.  20310 

Director,  Army  Ballistic  Research  Laboratory 
ATTN:  Mr.  J.  Meszaros 

Dr.  N.J.  Huffington,  Jr. 

Dr.  J.  Santiago 

Aberdeen  Proving  Ground,  MD  21005 

Commander,  U.  S.  Army  Missile  Command 
ATTN:  Dr.  R.  Rhodes 

Dr.  S.  Smith 
Huntsville,  AL  35809 

Commander,  Harry  Diamond  Laboratories 
ATTN:  AMXDO-RBF , Dr.  R.  Oswald 

AMXDO-NP,  Dr.  F.  Wimenitz 
2800  Powder  Mill  Road,  Adelphi,  MD  20783 
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Commander,  Picatinny  Arsenal 
1 ATTN:  Mr.  M.  Allen 

1 Mr.  M.  Weinstein 

1 Mr.  B.  Frank 

Dover,  NJ  07801 

Commander,  U.S.  Army  Combat  Development  Command 
1 ATTN:  Technical  Library 

Institute  of  Nuclear  Studies,  Fort  Bliss,  Texas  79936 

Commander,  Air  Force  Materials  Laboratory,  Air  Force  Systems  Command 
1 ATTN:  LNE/Major  H.  Keck 

1 LNC/Dr . D.  Schmidt 

Wright- Patterson  Air  Force  Base,  Ohio  45433 

Department  of  the  Navy,  Naval  Ordnance  Systems  Command 
1 ATTN:  ORD-03331 , Mr.  M.  Kinna 

Washington,  D.C.  20360 

Commander,  Naval  Surface  Weapons  Center 
1 ATTN:  Mr.  L.  Gowen 

1 Mr.  F.  Koubek 

Silver  Springs,  MD  20910 

Los  Alamos  Scientific  Laboratory 
1 ATTN:  GMX-6,  Dr.  J.W.  Taylor 

P.O.  Box  1663,  Los  Alamos,  NM  87544 

Space  and  Missile  Systems  Organization 
1 ATTN:  RSSE/Major  J.  McCormack 

P.O.  Box  92960,  World  Way  Postal  Center,  Los  Angeles,  CA  90009 

Sandia  Laboratories 
1 ATTN:  Dr.  D.  Munson 

1 Dr.  W.  Herrmann 

1 Dr.  L.D.  Bertholf 

1 Dr.  B.  Butcher 

1 Dr.  J.  Lipkin 

P.O.  Box  5800,  Albuquerque,  NM  87115 

Aerospace  Corporation 
1 ATTN:  Dr.  R.  Cooper 

1 Dr.  W.  Barry 

P.O.  Box  92957,  Los  Angeles,  CA  90009 

AVCO  Corporation,  Government  Products  Group 
1 ATTN:  Dr.  W.  Reinecke 

1 Mr.  P.  Rolincik 

201  Lowell  Street,  Wilmington,  MA  01997 

Bell  Telephone  Laboratories,  Inc. 

1 ATTN:  Mr.  E.G.  Denigris 

1 Mr.  M.F.  Stevens 

Murray  Hill,  NJ  07871 
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Effects  Technology,  Inc. 

ATTN:  Dr.  R.  Wengler 

Mr.  E.  Steele 

P.0.  Box  30400,  Santa  Barbara,  CA  93105 

Fiber  Materials,  Inc. 

ATTN:  Mr.  Maurice  Subilia,  Jr. 

Mr.  L.  Landers 

Biddeford  Industrial  Park,  Biddeford,  ME  04005 

General  Electric  Company,  Valley  Forge  Space  Technology  Center 
ATTN:  Mr.  K.  Hall 

Mr.  R.  Sullivan 
Mr.  J.  Braze  1 

P.0.  Box  8555,  Philadelphia,  PA  19101 

Kaman  Sciences  Corporation 
ATTN:  Mr.  F.  Shelton 

P.0.  Box  7463,  Colorado  Springs,  CO  80933 
Ktech 

ATTN:  Dr.  D.  Keller 

Mr.  N.H.  Froula 

911  Pennsylvania  Avenue,  N.E.,  Albuquerque,  NM  87110 

Lockheed  Missiles  and  Space  Company 

ATTN:  Mr.  D.  Aspinwall 

P.O.  Box  504,  Sunnyvale,  CA  94088 

Lawrence  Livermore  Laboratory 
ATTN:  Dr.  E.M.  Wu 

P.O.  Box  808,  L-421,  University  of  California,  Livermore,  CA  94550 


Martin  Marietta  Aerospace 
ATTN:  Mr.  M.  Hendricks 

Mr.  L.  Kinnaird 
Mr . F . Koo 

P.O.  Box  5837,  Orlando,  Florida  32805 


McDonnell  Douglas  Corporation 
ATTN:  Dr.  H.  Hurwicz 

5301  Bolsa  Avenue,  Huntington  Beach,  CA  92647 


Prototype  Development  Associates,  Inc. 

ATTN:  Dr.  J.I.  Slaughter 

Mr.  J.  Schutzler 

1740  Garry  Avenue,  Suite  201,  Santa  Ana,  CA  92705 


R&D  Associates 
ATTN:  Dr.  A.  Field 

525  Wilshire  Blvd.,  Santa  Monica,  CA 
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Southwest  Research  Institute 
ATTN:  Mr.  A.  Wenzel 

8500  Culebra  Road,  San  Antonio,  Texas  78206 


Stanford  Research  Institute 
ATTN:  Dr.  D.  Curran 

Dr.  L.  Seaman 

333  Ravenswood  Avenue,  Menlo  Park,  CA  90250 

Stanford  University,  Department  of  Applied  Mechanics 
ATTN:  Prof.  E.H.  Lee 

Stanford,  CA  94305 

TRW  Systems  Group 
ATTN:  Mr.  D.  Gamble 

One  Space  Park,  Redondo  Beach,  CA  90278 

Terra  Tek,  Inc. 

ATTN:  Dr.  A.H.  Jones 

420  Wakara  Way,  University  Research  Park,  Salt  Lake  City,  Utah 

Defense  Documentation  Center,  Cameron  Station,  Bldg.  5 
5010  Duke  Station,  Alexandria  VA  22314 

Director,  Army  Materials  and  Mechanics  Research  Center 


ATTN: 

DRXMR-H, 

Mr. 

J.  Dignam 

DRXMR-H, 

Mr. 

L.  Aronin 

DRXMR-H, 

Dr. 

S . C . Chou 

DRXMR-H, 

Dr. 

D.  Dandekar 

DRXMR-H,  Major  L.  Abramson 

DRXMR-AP 

DRXMR-PL 

DRXMR-PR 

Watertown,  MA  02172 
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